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1.  Introduction 


The  past  decade  witnessed  an  extensive  development  in  the  field  of  satellite 
radiotomography  (RT)  of  the  ionosphere,  which  significantly  extends  the 
capabilities  of  research  of  the  nearspace  environment.  Some  experimental  teams 
in  many  countries  have  been  working  on  the  problem  of  reconstruction  of 
ionospheric  sections  by  means  of  ray  RT  methods.  Since  1990  a  number  of 
experimental  studies  in  this  direction  have  been  reported  in  the  literature  [1-17]. 
Some  comparisons  with  the  data  from  incoherent  scattering  radars  and 
ionosondes  have  generally  shown  an  acceptable  quality  of  tomographic 
reconstructions  [12-14].  All  the  shown  experiments  were  realized  by  middle-orbit 
navigation  satellites  (  1000  km,  like  ’Transit"-USA  and  "Cikada''-Russia).  Surely 
the  development  of  these  data  by  means  of  high-orbital  navigation  satellites  is  of 
great  interest  (for  example,  like  "GPS"-USA  and  "GLONASS" -Russia  (  20,000  km)). 
The  consideration  of  a  combination  of  middle  and  high-orbital  satellites  is  usefull,  too. 
Such  radiotomographic  systems  allow  to  search  the  structure  of  magnetosphere, 
protonosphere,  etc,  and  to  study  the  influence  of  these  mediums  on  riavigation  and 
connection  systems  in  details. 

In  connection  with  the  experiments  performed  the  questions  arise  concerning 
the  data  quality,  the  accuracy  of  tomographic  reconstruction.  Given  strong  gradient 
of  electron  density,  one  of  substantial  contraints  imposed  on  ray  radiotomography  is 
the  ray  refraction.  Ignoration  of  refraction  limits  the  resolving  power  of  RT  systems. 
Taking  the  refraction  into  account,  it  allows  to  Increase  the  possible  resolving  power, 
but  it  transforms  to  more  complex  problems  of  nonlinear  RT.  Supposed  that  the  set  of 
integrals  over  retractable  rays,  which  curvature  is  determined  by  the  propagating 
medium  are  known.  That  is  why  the  exploration  of  abilities  of  nonlinear  RT  is  seems 
important. 

The  aim  of  this  work  is  to  develop  algorithms  and  programs  for  the  ray  nonlinear 
(with  refraction)  Radiotomography  of  the  Ionosphere  and  Radiotomography  using 
middle  and  high-orbital  satellites. 
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Algorithms  and  Programs  for  Radiotomography  using  middle  and  high- 

orbital  satellites. 


1.1.  Description  variants  of  sateliite  Radiotomography  using  niiddie  and  high- 

orbital  satellites  of  opportunity. 

As  modern  navigation  systems  on  high-orbital  satellites  (  like  "GPS"  -  USA  and 
"GLONASS"  -  Russia)  are  operated  on  a  frequency  of  1  GHz,  in  this  case  it  is  enough 
to  consider  problems  of  linear  tomography  ignoring  the  ray  refraction.  The  ray 
refraction  is  a  factor  at  a  frequency  of  hundreds  of  MHz.  Such  problems  of 
radiotomography  would  be  considered  in  the  second  part.  Modern  navigation  systems 
on  high-orbital  satellites  mentioned  in  the  first  part  (  "GPS"  or  "GLONASS"  )  are 
supposed  to  operate  on  frequences  nearly  1 .2  and  1 .5  GHz. 

Firstly  let  us  remind  the  main  relationships  which  connect  the  experimentally 
measured  physical  values  with  the  parameters  of  propagating  medium.  Known  is  [18- 
20]  the  differential  phase  between  the  two  L-band  carriers,  at  /j  =  1 .575  GHz  and 
/2  =  1 .227  GHz  can  be  used  to  make  precise  relative  total  electron  content  (TEC) 
measurements.  The  differential  group  delay  between  the  10.23  MHz  modulation  on 
the  two  L-band  carriers  is  used  to  fit  these  precise  relative  TEC  measurements  to 
an  absolute  scale. 

The  difference  of  simultaneous  dual-frequency  GPS  code  measurements  Pj  , 
Pj  (in  meter)  equals  the  TEC  along  the  transmitter-receiver  ray  [18-20]. 

TEC  s  J  Nda  =  -  Pi)  +  (1 ) 

Where  J  Nda  denotes  the  integration  of  the  electron  density  N  along  the  signal 
path.  Merging  all  bias  terms  (differential  instrumental  group  delays  in  the  receiver 
and  in  the  satellite,  and  multipath  and  observation  noise)  yields  the  error 

2 

The  scaling  factor  S  converts  units  of  distance  [m]  to  units  of  TEC  (electrons/?W  ) 
and  is  derived  from  the  GPS  frequencies  with: 


s 


s  = 


0.025 


9.52*10‘®m'^ 


The  observation  equation  for  the  TEC  from  phase  measurements  (phase  puth) 
OjjOj  [m]  iooks  similar  [18-20].  The  error  include  the  differential 

instrumental  biases  of  the  phases  and  additional  bias  term  exist  due  to  the  carrier 
phase  ambiguities 


TEC  =  j  Nda  -  S{^2  ~  ‘^i)  +  -  M2A2)  +  (2) 

Where  the  \,A2  denotes  the  wavelengths  of  sounding  waves. 

Thus  phase  and  group  delay  measurements  by  means  of  modern  navigation 
systems  GPS  and  GLONASS  allow  to  find  linear  integrals  of  the  electron 
concentration.  The  following  fact  is  of  great  importance  for  further  consideration  of 
radiotomographic  problems  by  means  of  GPS  and  GLONASS.  TEC  can  be  measured 
to  a  relative  accuracy  of  approximately  3*10^'*  e/  /  m^,  and  absolute  values  can 
be  determined  to  approximately  (l-2)*10^®e// m^,  plus  any  unknown  GPS 
satellite  differential  code  offsets  [21-24].  In  other  words  given  methods  allow  to 
determine  relative  TEC  precisely  enough,  but  they  showed  unacceptable  mistake  in 
determination  of  absolute  TEC.  This  case  is  similar  to  the  situation  in  Transit-Cikada 
RT  systems  where  major  difficulty  is  that  the  linear  integral  of  the  electron 
concentration  is  proportional  to  the  absolute  phase  3>,  whose  accurate  determination 
is  practically  impossible,  but  the  differences  of  linear  integrals  (  phase  difference  ) 
are  measured  with  high  precision.  That  is  why  the  new  approach  was  taken  earlier  to 
determine  the  difference  of  linear  integrals  [1-3,  25],  which  (approach)  brought 
satisfyed  results  [1-3,  13,  14,  17].  Further  the  same  approach  would  be  taken  to 
determine  the  difference  of  linear  integrals. 

Let  us  consider  possible  variants  of  using  of  midlle  and  high-orbital  satelittes 
for  radiotomographic  purposes.  Mostly  geometry  of  experiments  is  essential  for 
qualitative  radiotomographic  reconstruction:  the  correlation  of  receivers  and 
transmitters,  the  way  of  crossing  the  reconstrtuctive  area  by  rays,  angle  between  rays. 
The  program  "SATTOMO"  was  worked  out  to  make  geometrical  analysis  of 
registration  scheme.  Given  program  allows  to  analyze  radiotomographic  registration 
schemes,  which  use  ground  receivers  and  measuring  systems  on  satelittes  with  cyclic 


orbits,  such  as  Transit,  Cikada,  GPS,  GLONASS,  Microlab,  etc.  The  description  of  this 
program  is  cited  in  the  end  of  this  article. 

Heights  of  orbits  of  sateiiites  in  program  SATTOMO  can  be  varyed,  respectively 
the  period  and  the  rotating  angular  velocity  is  changing.  Zero  height  is  corresponded 
to  receivers  on  the  Earth,  whose  angular  velocity  is  supposed  to  be  zero.  Given  critical 
altitude  is  the  minimum  height  of  passage  of  rays  over  the  Earth  .  This  height  can  be 
chosen  equal  to  zero.  In  other  cases  the  choise  of  this  height  limits  the  range  of  ray 
registration,  and,  for  example,  limits  the  rays  passing  through  lower  ionosphere.  If  the 
number  of  satellites  is  given,  it  is  possible  to  determine  the  direction  of  satellites 
rotation  in  any  of  two  orbits  (+  counterclockwise,  -  clockwise)  Angular  location  in  the 
initial  moment  of  time  is  given  in  degrees  in  polar  coordinate  system  where  the  zero 
angle  is  corresponded  to  horisontal  X-coordinate.  Program  also  includes  the  setable 
time  of  tomographic  measurments  in  seconds,  time  interval  between  drawing  of 
two  neiboring  rays  (convenient  time  is  30-60  s  ,  if  the  less  time  interval  is  chosen,  the 
picture  of  rays  will  be  too  densely)  and  time  of  drawing  of  the  total  picture  on  the 
screen  in  seconds.  The  last  is  wished  to  be  little  (5-10  seconds)  to  have  an 
opportunity  to  examine  the  satellite  motion. 

The  use  of  the  only  high-orbital  navigation  system  like  GPS-GLONASS  with  the 
ground  receivers  for  purposes  of  radiotomography  of  ionosphere  and  protonosphere 
appears  to  be  unperspective  with  the  respect  to  little  angular  velocity  of  high-orbital 
satellites.  Indeed,  for  example,  the  GPS  angular  velocity  is  equal  to 
CO  =  0.0001456  rad  /  sec,  what  is  seven  times  less  than  the  middle-orbital  satellite 
angular  velocity  equal  to  0.0009966  rad  /  sec.  Such  satellites  pass  small  distances  in 
tens  of  minutes  (20-40)  and  fill  a  small  angle.  The  use  of  a  greater  time  intervals  is 
almost  impossible  because  of  the  time  change  in  ionosphere.  In  fig. la  shown  -  is 
the  angle  of  filling  the  space  by  rays  over  the  ground  receiver  during  the  GPS  fly  over 
the  head.  This  angle  (A0  =  Aa(l-i-i?/  H))  is  close  to  polar  angle  of  satellite  Aor, 
as  the  height  of  satellite  is  H  =  20200  km  what  is  essentially  exceeds  the  Earth  radius 
R  =  6370  km.  Polar  angle  is  changes  proportionally  to  the  time  Acr  -  cot.  For 
example,  in  one  hour  (  equal  to  1/12  of  period  of  rotating)  the  polar  angle  will  change 
only  for  2n  j  Yl-  30°.  The  intersection  of  the  rays  (fig.  1b)  for  receivers  located  at  a 
distance  of  Ar  (  on  the  Earth)  occurs  on  the  height 

hQ'HiAT / 2tan(A6 / 2)f>iAT / 2tan(cot / 2) ,  what  is  (  for  t^lSOOs)  «  12Ar  /  n.  If 

the  distance  between  receivers  is  greater  than  250  km,  than  the  rays  of  neiboring 
receivers  are  intersected  on  heights  more  than  1000  km.  To  make  the  rays  intersect 
in  the  lower  boundary  of  ionosphere  (150  km),  receivers  should  be  located  densely 
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(At  <  60km).  Fig. 2  describes  the  geometry  of  radiotomographic  experiment  with 
one  and  four  receivers  on  the  Earth  (  angular  location  30°  ,40°  ,50°  ,60°  ).  Even  after  an 
hour  of  measurings  the  essential  "holes"  were  taken  place,  where  the  rays  didn’t  get 
into.  The  chosen  distance  between  receivers  was  too  great  to  examine  the  details  in 
such  a  scale.  Picture  didn't  change  high-quality  ,  as  the  receivers  have  been  coming 
together,  only  the  size  and  height  of  "holes"  were  reduced.  Geometry  of  experiments 
appears  to  be  unsuitable  for  qualitative  radiotomographic  reconstruction  with  the  use 
of  high-orbital  satellites  and  receivers  on  the  Earth.  Next  two  pictures  show  bad 

geometry  and  the  existence  of  "holes"  for  four  ground  receivers  (  7^  =3600  s  )  and 
for  two  GPS  satellites  (fig.3).  The  same  for  two  GLONASS  satellites  (fig.4).  However 
the  use  of  network  of  ground  receivers  on  the  given  square,  each  of  that  registrates 
signals  simultaniously  from  GPS  and  GLONASS  satellites  in  various  directions  is  waited 
to  be  perspective.  Such  a  system  allows  to  determine  parameters  of  ionosphere  in 
three-dimensional  area  over  the  given  region. 

The  using  of  combination  of  different  satellite  systems  appears  to  be  more 
perspective  for  radiotomography.  Similar  satellite  systems  were  already  worked  out  for 
scientific  investigations.  Recently,  on  April,  3,  1995  Orbital  Sciences  Corporation 
successfully  launched  the  Microlab- 1  scientific  spacecraft  (with  altitude  600  km). 
However  the  only  middle-orbital  satellite  (  with  altitude  equal  to  500-2000  km  )  also 
does  not  give  the  good  geometry  of  experiment.  Fig. 5-6  show  the  schemes  of 
experiments  with  one  middle-orbital  satellite  (T^  =  I8OO5)  in  combination  with  high- 
orbital  satellites:  one  Microlab  with  two  GPS  (angular  location  of  GPS  in  the  initial  time 
moment  -  20,  110  degree)  and  Microlab  (  25  degree)  (pic.5),  and  one  Transit  (  30 
degree)  with  three  GLONASS  (angular  locaton  in  the  initial  time  moment  -  15  ,  60  , 
105  degree)  (pic. 6).  These  pictures  show  slight  intersection  of  rays.  With  the  help  of 
such  a  scheme  it  is  possible  to  explore  only  part  of  protonosphere,  where  rays  are 
intersecting,  under  the  condition  that  this  area  has  strong  enough  local  disturbance 
and  the  contribution  of  inside  areas  of  protonosphere  is  negligible. 

The  geometry  of  experiment  is  essentially  better  if  the  combination  of  a  few 
middle-orbital  satellites  on  the  orbit  with  high-orbital  satellites  is  used.  Fig.7  gives  the 
scheme  of  an  experiment  with  three  satellites  Transit  (  angular  location  of  the  initial 
time  moment  -  0,  30,  60  degree)  and  one  GPS  (15  degree).  It  is  seen  that  thirty 
minutes  measurements  allow  to  reconstruct  the  structure  of  protonosphere  on  a  wide 
enough  band.  Fig.8  describes  the  scheme  of  radiotomographic  registration  with  three 
Transit  (angular  locations  in  the  initial  time  moment  0,  30,  60  degree)  and  three 
GLONASS  satellites  (angular  locations  in  the  initial  time  moment  15,  60,  105  degree). 
Here  thirty  minutes  measurements  allow  to  reconstruct  the  structure  of  protonosphere 
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in  three  intersective  bands.  Though  this  geometry  is  worse  than  a  typical  geometry  of 
experiments  on  satellites  Cikada  and  Transit  [1-17],  as  the  range  of  angles  of 
intersection  between  rays  is  smaller,  nevertheless,  as  it  would  be  shown  lower,  it 
allows  to  get  radiotomographic  cross-sections  of  protonosphere. 

It  is  a  remarkable  fact  that  the  use  of  several  middle-orbital  satellites,  rotating  in 
the  opposite  direction,  gives  a  good  geometry  with  a  set  of  fan  beams.  Fig.9  shows 
the  scheme  of  experiment  with  one  GPS  and  three  Microlab  satellites  (angular  location 
in  the  initial  time  moment  90,  120,  150  degree),  rotating  in  the  opposite  to  GPS 
(80  degree)  direction.  After  the  registration  period  of  twenty  minutes  =  12005)  the 
ray  structure,  which  consists  of  three  beams  is  seen  well.  Fig.  10  shows  one  of  such  a 
beam  with  the  height  of  "focus"  multiplicity  of  12000  km  for  one  Microlab  satellite.  Let 
us  cite  one  more  interesting  example  of  using  of  combination  of  middle  and  high- 
orbital  satellites.  Fig. 11  shows  the  scheme  of  fifteen  minutes  (T^  =  9005)  measurings 
between  one  Transit  satellite  (angular  location  in  the  initial  time  moment  90  degree) 
and  one  GPS  (25  degree).  Such  satellite  orientation  allows  to  fill  the  area  of  nearearth 
ionosphere  by  quasiparallel  rays.  Fig.  12  describes  analogious  scheme  (T^  =  9605) 
with  one  Transit  satellite  (  angular  location  in  the  initial  time  moment  90  degree)  and 
two  GLONASS  (0,  45  degree).  This  geometry  allows  to  fill  the  area  around  all  the 
Earth  by  quasiparallel  rays. 

First  of  all  we  want  to  remind  the  geometrical  equations,  which  connect  the 
transmitters  and  recievers  locations  [3,6,8].  Here  we  will  conditionally  indicate  the 
locations  of  transmitters  and  recievers  on  the  concrete  orbits,  meaning  that  the 
observing  result  will  not  depend  on  changing  this  locations.  We  introduce  a  series  of 
parameters  characterizing  the  geometry  of  the  recording  scheme  in  polar  coordinates 
(r,  a).  In  Fig.  13  -  are  the  coordinates  of  one  of  the  receivers  located  on 

the  first  orbit  r  =  (also  the  reciever  can  locate  on  the  Earth  surface  -  R^  =  R  - 
Earth  radius);  (i?2,ao)  '  ^^e  the  coordinates  of  the  satellite  with  transmitter  on  the 
second  orbit  /*  =  elevation  of  the  i?2  satellite;  ^  =  -  n  !  2 

is  the  angle  to  the  satellite  measured  from  vertical;  O  is  the  center  of  the  earth;  axis 
0-0'  is  axis  of  the  polar  coordinate  system.  We  consider  here,  that  satellites  rotate  in 
single  or  nearer  planes,  what  is  important  for  RT  purposes.  On  the  basis  of  simple 
geometrical  relations  for  any  point  in  the  ionosphere  with  coordinates  (r,  a )  located 
at  the  distance  of  I  from  the  receiver  the  following  equations  are  satisfied: 


I  r  i?i 

sin(a,.  -  a)  sin(;r  f  2  +  p)  sin(;r  -{n  /  2  +  p)-ia^  -  a)) 

From  this,  we  obtain  an  equation  for  r{a)  of  the  straight  ( =const)  ray 

r{a)  =  i?i  cos(;0 )  /  (cos(^  +  a^  -  a)) 

and  the  relation  which  is  inverse  of  it 


a{r)  =  a^  +  p  -  arccos(-^  -  cos  p ) 


The  relation  between  p  and  a ,  r  follows  from  (4); 


tan  p  =  (cos(a;  -  a)-  /  r)  /  (sin(a,  -  a)) 


Using  (5),  we  obtain  a  formula  for  an  element  of  the  ray  of  length  da: 


da^  =  [I +  r\—y}^dr^ 


dr'  r^  -  Rl  cos^  p 


Then,  the  relation  for  the  measured  linear  integrals  (1,2)  with  respect  to  the 
electron  concentration  wili  have  the  form 


f  N{r,a)rdr 

H  d2 


r  -  R{  sin  q) 


In  the  place  of  the  polar  coordinates  (r,  a),  hereafter  it  is  convenient  to  use 
the  orthogonal  system  {h,T):h  =  (r  -  R^)  is  the  height  above  the  first  orbit  or  above 
the  earth's  surface  and  t  =  aR^  is  the  "transverse"  (horizontal)  distance  along  first 
orbit.  Here,  ray  equation  (4)  wili  no  longer  be  a  straight  line 


COsB 

h{T)  =  R[ - - 1] 

1  cos{p  +(j^-t)/ Ri) 


Here,  (r,,/i  =  0)  are  the  coordinates  of  the  receivers.  The  relation  which  is  the 
inverse  of  (8)  is  similar  to  (5). 


r  -  r,  =  -  arccos(-^^cos^ )] 


(9) 


In  this  case,  subject  to  (7),  the  linear  integrals  of  type  (1-2)  have  the  form 


r/iQ  FQi,  r){R,+  h)dh 


(10) 


Integration  with  respect  to  the  ray  connecting  receiver  i  (ti  =  CCjRj)  to  the 
satellite  is  replaced  according  to  (7)  by  integration  with  respect  to  the  height  above 
the  first  orbit  to  the  height  of  the  satellite  .  The  elevation  p  is  determined  (6)  by 
position  of  the  satellite  (/io,7o).  The  linear  integral  I{p,Ti)  is  dependent  on  the 
coordinate  of  the  receiver  r,-  and  the  elevation  of  the  satellite  p  .  The  linear  integral 
is  the  TEC  here,  the  reconstructed  function  F  will  be  proportional  to  N.  Since  there 
cannot  be  a  large  number  of  receivers  and  the  range  of  angles  p  is  limited,  it 
is  inadvisable  to  examine  methods  for  analytical  inversion  of  such  linear  integrals 
and  methods  for  integral  transformations.  In  given  case  small-positions 
(smallforeshortening)  tomography  is  intended  from  the  start  to  solve  the  problem  in 
discrete  form  and  to  use  algebraic  reconstruction  algorithms  or  methods  for 
expansion  into  finite  series. 


Description  of  the  the  programm  “sattomo”  for  vizuaiisation  geometry  of 

satellite  RT  experiments. 

This  program  allow  to  view  on  the  screen  the  different  versions  of  the 
radiotomography  measurements  geometry  for  the  one-  and  two-orbit  satellite 
configurations. 

Input  parameters  and  files. 

Command  line  format  is:  SAT2.exe  <config_file_name>  [p].  First  parameter  -  the 
name  of  file  that  defines  all  the  input  parameters  for  RT  system  geometry. 


Second  parameter  (optional),  if  set,  allow  to  print  geomerty  picture  on  the 
EPSON-compatible  printer. 

The  configuration  file  has  the  following  format:  several  lines  started  at  1st 
position  in  form:  PAR1,PAR2,....  [;  comments]: 

Line  1:  HJow  H_high  H_crit 
Line  2:  NSatl  [+  or  -]  NSat2[+  or  -] 

Lines:  Pos1(1)  Pos1{2) .  Pos1{NSat1) 

Line  4:  Pos2(1)  Pos2(2) .  Pos2{NSat2) 

Line  5:  PassTime  Step  DrawTime 

The  above  mentioned  parameters  are: 

-  HJow,  H_high  -  orbit  heights  (kilometers)  of  2  groups  of  satellites.  HJow=0 
defines  the  receiver  on  the  Earth's  surface.  -  H_crit  -  minimum  height  for  the  ray 
traectory.  If  H_crit=0,  the  ray  may  touch  the  Earth's  surface. 

-  NSatl ,  NSat2  -  the  number  of  satellites  on  the  heights  HJow,  H_high.  The  signs 
'+'  or  '-'  are  optional,  but  if  they  are  set,  must  follow  immediately  after  this 
numbers  and  define  the  direction  of  the  satellite  moving  (clockwise  or  opposite). 
If  no  sign  character  set,  '+'  assumed. 

-  The  arrays  of  Post  and  Pos2  define  the  initial  positions  of  satellites.  Must  be  set 
in  degrees. 

-  PassTime  -  full  time  of  real  measurements. 

-  Step  -  the  ray  drawing  step  time. 

-  DrawTime  -  time  that  needs  to  draw  picture  on  the  screen. 

The  tree  last  parameters  must  be  set  in  seconds. 


1.2.  The  solution  of  the  direct  problem  of  TEC  calculation  for  high-orbital 

satellites  Radiotomography. 

The  purpose  of  this  section  is  to  describe  the  technique  and  program  for 
solving  the  direct  problem  of  ionosphere  radio  probing,  i.e.  the  problem  of 
obtaining  the  TEC  and  the  TEC  difference  from  the  data  on  the  electron  density 
distribution.  In  view  of  calculation  the  direct  problem  solution  amounts  to 
calculatiry  the  integral  (10),  i.e.  the  TEC  or  the  TEC  difference  of  such  integrals  for 
a  small  variations  of  the  satellite  positions.  For  the  computer  modeling  of  RT 


problems  it  is  necessary  to  calculate  a  series  of  such  integrals  for  arbitrary 
positions  of  the  receivers  and  the  transmitters  on  the  satellites,  therefore  we 
accomplished  the  program  for  calculating  TEC  and  the  TEC  difference  for  arbitrary 
positions  of  the  receivers  and  the  transmitters  on  the  satellites. 

Since  the  integrand  has  no  peculiarities,  the  calculation  can  be  performed 
using  the  rectangle  technique  or  the  Simpson  method.  Let  us  evaluate  the  necessary 
step  of  numerical  integration  and  the  accuracies  obtained  in  this  way.  it  is  well 
known  that  errors  of  numerical  integration  by  the  rectangle  method  and  the 
Simpson  method  are  equal  to: 

^3  pm  ^5  pm 

"  I2m^  ’  "  2880m^ 

where  S  is  the  integration  interval  length  (the  distance  between  satellites),  m  is  the 
number  of  the  integral  discretization  elements.  Integration  errors  made  using  the 
rectangle  and  Simpson  methods  are  proportional  to  the  values  of  the  second  and 
fourth  derivatives,  respectively,  of  the  integrand  at  a  certain  point  within  the 
integration  interval.  For  aii  approximate  estimation  of  errors  it  is  sufficient  to 
evaluate  the  second  and  fourth  derivatives  as  a  result  of  dividing  a  characteristic 
value  of  the  function  F  by  the  square  or  the  fourth  power  "a",  respectively,  of  the 
structural  irregularities  of  F.  Limitations  associated  with  diffraction  effects  as  well  as 
those  of  the  linear  tomography  problem  make  it  impossible  to  reconstruct  details 
smaller  than  10-20  km  using  the  method  of  linear  ray  RT  [3,8],  therefore 
fl  >  10  -  20km  The  value  of  m  is  equal  to  the  result  of  dividing  S  by  the 
integration  step  As.  Hense  the  estimations  of  absolute  and  relative  errors  are: 

SF  SF 

(12) 

^^^s(AA/«)V2880. 

When  the  integration  step  As  =  0.5km  and  a=20km  the  relative  error  of  the  rectangle 
method  is  0.5*10'“’  (for  the  Simpson  method  it  is  smaller  than  10'”)  which  is 
quite  satisfactory  for  RT  applications. 


Description  of  the  program  for  calcuiation  TEC  and  TEC-difference  for  the 

high-orbit  RT  measurements 


System  Requirements 


-  Computer:  IBM  AT-486  or  compatible  (with  coprocessor) 

-  Operating  System:  MS-DOS  or  PC-DOS  version  3.0  and  later 

-  Memory:  at  least  Extended  memory  8  Mbytes 

(depends  on  geometry  and  type  of 
approximation  reconstructed  function) 

-  Hard  Disk  Space:  35  Mbytes 

-  Software:  NDP-FORTRAN-486  Compiler  Ver.3. 1 

NDPLink  Ver.3.0  (C)  MicroWay,  Inc. 


1.  Program  <dir_rth1.for> 

This  program  solves  direct  problem  for  high-orbit  RT  measurements, 
namely,  determines  the  model  structure  and  calculates  the  TEC  or 
TEC-difference  integrals  on  the  model  structure. 

Input  parameters  and  files: 

When  starting  program,  you  may  choose  the  type  of  RT  measurements 
(TEC  or  TEC-difference).  Then  program  asks  about  inegrals  to  calculate. 

If  you  need  only  model  structure  you  can  answer  'No'. 

Input  file  <task.dir>  contains  information  about  RT  measurements 
geometry  and  defines  the  names  for  output  files.  The  example  of  this 
file  is  enclosed.  It  contains  the  following  lines: 

'model. fun'  -  input  file  that  describes  model  structure 
'model.grd'  -  output  file  with  model  structure  (GRD-format) 

NF  NR  -  number  of  discrete  on  the  horizontal  and  vertical  grid 
0.  H_GPS  -  high  orbit  altitude  (GPS  satellites)  over  low  orbit  in  km 
TFIST  TFIN  -  size  of  reconstructed  area  in  km  on  the  low  orbit  height 
STEPS  -  number  of  the  discrete  steps  along  the  ray 
Nrec  -  number  of  receivers 

LOs  LOe  HOs  HOe  Nrays  -  Nrec  lines  with  initial  and  final  positions  of 
satellites  on  the  low  and  high  orbits  in  km  and  number  of 


rays  for  this  couple  of  satellites 
’linjnt'  -  Nrec  lines  with  the  names  of  output  files  for  integrals 
for  each  receiver 

The  model  structure  is  determined  by  function  <FMODEL>  which  uses 
functions  <HOMCOS>.  Function  parameters  are  defined  in  file  <model.fun>. 
The  example  of  this  file  is  enclosed. 

Output  files: 

MODEL.GRD  -  file  with  model  structure  (GRD-format) 

Files  <linjnt>  -  arrays  of  either  TEC  or  TEC-difference  integrals 

Compilation 

mf77  dir_rth1.for  -ol  -486 
RUN 

ndprun  dir_rth1  .Iti 


1.3.  The  construction  of  the  projection  operators  for  the  high-orbital  RT 
measurements  using  various  approximation  methods. 

In  this  section  we  consider  the  various  methods  to  constuct  the  projection 
operator  L  that  translates  the  sought  function  of  the  electron  concentration 
distribution  F  into  a  set  of  integrals  (TEC)  I: 


LF  =  I 


First  of  all  we  shell  perform  the  discretization  procedure  for  equations  (1)  to 
prepare  the  numerical  calculations.  We  shall  use  the  orthogonal  coordinate  system 
where  h  -  is  the  height  above  the  low  satellite  orbit  and  t  -  is  transverse 
distance  along  this  orbit.  To  digitize  the  sought  function  F(h,r)  in  a  fixed 
rectangular  grid  we  divide  the  rectangular  reconstruction  region  into 

/Mq  heights  (m  <  a/Iq  )  and  horizontal  samples  (1  <  n  <  «o  )  and  replace  function  by 
a  piecewise-constant  approximation,  or  to  represent  F  by  a  system  of  (mo*«o)  basis 


functions  equal  to  unity  in  certain  rectangle  and  zero  in  all  others.  The  value  of  the 
function  F{h,T)  in  a  fixed  (//?''«)  rectangle  we  define  as 

We  perform  digitization  of  the  linear  TEC  integrals  (10)  according  to 

the  initial  direction  to  the  satellite  on  the  high  orbit  from  receiver  number  i 
placed  on  the  low  orbit  satellite  with  current  coordinate  Tj  .  The  set  of  elevation 
(6)  of  all  moving  receivers  define  a  series  of  discrete  values  of  the  linear  integrals 
(10)  ly  =  I(Py,Ti),  where  i  -  receiver  number,  j  -  ray  number  from  set  for  receiver  i. 

The  number  of  rays  is  determined  by  the  parameters  of  the  recording  system.  The 
point  in  the  rectangles  at  which  the  samples  of  F{h,  t)  are  selected  is  not  especially 
important;  this  may  be  at  the  middle  of  the  rectangles  or  nodes  of  the  grid. 

The  simplest  method  of  project  operator  discretization  is  piecewise-constant 
approximation,  when  the  coefficients  of  matrix  A  are  proportional  to  the  ray  length  in 
elements  of  the  discrete  grid.  Designating  the  length  of  ray  (i,])  in  cell  (m,n)  as 
,  we  obtain  the  system  of  linear  equations 


or  after  "renumbering”  of  the  ray  (i,j)  -  J  and  the  cells  of  the  ionosphere  (m,n)  -  M 

LfF^=Ij  (13) 

System  (13)  may  be  either  overdetermined  or  sub-definite.  The  problem  of  high- 
orbit  tomographic  reconstruction  according  to  TEC  measurements  is  to  determine  the 
set  of  discrete  samples  {F^}  in  the  known  grid  according  to  the  set  {ly}- 

This  simple  method  of  projection  operator  constructing,  however,  gives  the 
high  error  value  of  the  direct  problem  solution  (particulary  for  the  small  number  of 
discrete  elements).  Besides  that  the  piecewise-constant  approximation  is  invalid 
for  the  problem  of  high-orbit  RT  according  to  TEC-difference  measurements.  The  fact 
is  that  the  data  here  will  be  derivatives  of  linear  integrals  of  type  (10): 
D  =  dl  /  da^,  or  finite-difference  ratios  of  the  increment  A/  of  the  linear  integrals  to 
the  increment  Aa,,  of  the  satellite  coordinate.  The  difference  of  TEC  measured  in 
the  experiment  can  be  determined  by  the  TEC  (1)  derivative.  Because  the  satellites 
move  uniformly  along  a  circular  orbits  it  is  possible  to  express  the  difference 
of  TEC  by  means  of  the  derivative  with  respect  to  the  angle  of  the  satellite  on  the 


high  orbit  ag  and,  hence,  TEC-difference  tomography  data  are  proportional  to 
A//  AcTo-  The  derivatives  of  the  linear  integrals  in  a  piecewise-constant 
approximation  of  the  sought  function  F  will  be  discontinuous.  This  is  because  each 
linear  integral  Is  the  sum  of  integrals  over  the  set  of  cells.  As  the  elevation  of  the 
high  orbit  satellite  changes,  the  ray  encounters  a  new  cell;  the  integral  with 
respect  to  every  of  this  cells  is  a  continuous  function  of  the  angle  of  the  satellite 
Uq  ,  but  the  derivative  of  the  linear  integral  with  respect  to  will  contain  a 
discontinuity  when  the  ray  contacts  the  corner  of  each  cell.  Therefore,  the 
piecewise-constant  representation  of  the  function  to  be  reconstructed  is  not 
appropriate  to  solve  the  TEC-difference  problem  correctly. 

To  ensure  continuity  of  linear  TEC  integrals  with  respect  to  the  coordinates  of 
the  both  satellites  (or  elevation  P  )  the  matrix  Ljj^  for  transition  from  the  function 
to  be  reconstructed  to  linear  integrals  should  be  calculated  differently.  The  main  idea 
of  other  methods  of  RT  operators  (matrices)  design  is  the  increase  of  the 
approximation  order  for  the  reconstructed  function  and,  hence,  for  the  matrix  L. 
We  used  the  piecewise-planar  approximation  based  on  the  triangular  elements, 
product  of  linear  approximations,  product  of  local  spline  approximations  and  the 
modification  of  the  last  method  also  including  derivatives  of  the  function  on  the 
discrete  grid.  The  accuracy  of  the  direct  problem  solution  is  increased  with  the 
higher  approximations  and  the  distinction  between  constructed  and  the  real 
"natural"  operators  is  reduced.  The  matrix  of  the  direct  problem  constructed  by 
these  methods  Ij  is  continuous  with  respect  to  the  angle  of  the 

satellite  on  the  high  orbit  ,  hence,  in  place  system  (13)  it  is  possible  to  obtain  a 
system  for  TEC-difference  data  by  differentiatiiig  (13)  with  respect  to  the  angle  cTo: 

Here,  Dj  =  AIj  /  Aa^  are  TEC-difference  data  and  =  ALj^  /  AcTo  is  finite  - 
difference  ratio  (or  derivative)  of  the  matrix  Ljj^  to  the  increment  of  the  angle. 

In  this  section  we  consider  examples  of  constructing  smooth  projection 
operators  of  the  direct  problem  wich  are  smooth  by  the  satellite  angle.  The  examples 
of  constructing  the  Ljj^  matrices  of  the  transition  from  the  reconstructed  function 
to  the  linear  integrals  (matrices  of  projection  operators)  are  useful  for  both  TEC- 
difference  RT  and  for  the  TEC  RT  but  the  first  is  preferable.  The  advantages  of  the 


first  method  was  repeatedly  shown  on  the  example  of  phase-difference  RT 
[6,8,25]. 

One  must  construct  such  a  Ljj^  matrix  that  would  provide  the  smooth  of 
linear  integral  over  the  satellite  angle  using  the  piecewise-planar  approximation  of 
the  sought  function  on  the  basis  of  triangular  elements.  Such  approximation  ensures 
the  smoothness  of  the  matrix  Lj^f  to  be  constructed.  The  smooth  function  FQi^r)  is 
replaced  by  a  continuous  polyhedral  approximation  surface,  according  to  which 
the  derivative  with  respect  to  the  satellite  angle  of  the  linear  integrals  is  already  a 
continuous  function.  Triangular  elements  are  obtained  naturally  from  a  grid  of 
rectangles:  each  cell  (m,n) 


=  A7-*A/i 

is  divided  by  a  diagonal  running  downward  and  left  -to-right, 

into  two  triangular  elements:  the  "lower"  and  "upper"  elements.  The  function 
F(h,T)  within  each  triangular  element  is  replaced  by  linear  approximation 

F(h,T)  =  a  +  br  +ch  (15) 

The  values  of  the  coefficients  (a,  b,  c)  in  each  finite  element  are  determined  from 
system  of  three  equations  for  three  boundary  points.  It  is  simple  to  immediately 
write  expressions  for  the  function  presentation  in  the  lower  (m,n)  element 

F(h,T)  =  F„  +  (h-h.)  (16) 


and  in  the  upper  (m,n)  element 


Fih,T) 


Ar 


+ ^ (17) 


As  before,  to  simplify  the  notation  we  will  renumber  the  values  of  the  samples  below: 

(m  +  l,n)  (M  +  l),(m,n  +  l)  (M  +  AM),(m  +  l,n  +  l)  (M  +  AM  +  1), 
where  AM  is  the  number  of  cells  horizontally  in  one  row.  Now  we  can  define  the 
result  of  integration  of  such  an  approximation  in  lower  element  M 


J  w(h)Fdh  -  JqFj^  +  /j.  ) 


(18) 


and  in  upper  triangular  element  M 

J  w(h)Fdh  =  +  J ^  (^Af+AM+l  “  '^Af+AAf  )  +  •^a(^m+am+i  “  -^M+i)  (^^) 

Here 

^  \y>mm  -  r.lrfA,  /;  =  ^JwWIrW  -  r..,]dA, 

(20) 

A  =^l«'W|A-A.lrfA./»  rw(A)[A-A.„ldA 

where  wifi)  =  (i?  +  sin^  F  +  2i?/i  +  Fih,r)  is  represented  in  the  form 

of  piecewise -planar  approximations  (16), (17)  in  each  finite  element,  /i„  -  lower 
boundary  of  the  cell  in  row  n,  h  -  height,where  the  ray  leaves  the  lower  element  and 
enters  into  upper.  After  integration  (18)  with  respect  to  ray  J  in  the  lower  element 
M  the  value  ~^h)  '®  entered  into  the  coefficient  Ljj^  ,  since  it  is  a 

coefficient  for  Ff^.  Correspondingly,  is  entered  into  and  into 

The  complete  value  of  the  coefficient  is  the  sum  of  integration  result  with  respect 
to  the  ray  J  in  six  around  triangular  elements.  The  integrals  with  respect  to  all  rays 
of  type  (20)  can  be  calculated  by  various  numerical  methods,  and,  in  each 
integration  step  idt  it  is  necessary  to  verify  that  the  ray  does  not  exceed  the  limits  of 
the  finite  element.  After  complete  numerical  integration  with  respect  to  all  rays, 
we  obtain  the  matrix  Lj^^.  The  matrix  Lj^f  is  related  to  the  set  {cTo}  of  positions  of 
the  satellite  on  the  high  orbit  and  the  corresponding  series  of  rays.  To  determine  the 
matrix  for  TEC-difference  tomography  problem  Aj^f  one  must  calculate  the  matrix  L' 
for  another  set  of  close  positions  of  the  satellites  with  a  fixed  increment 
{oTo  +  Atto)  for  the  first  satellite  and  {a^  +  Aaj)  for  the  second  one  and  to  find  the 
difference  between  them:  Ajj^f  =  {L /  Aoto-  Because  all  circular  satellites 
positions  are  proportional  to  time  it  is  equivalent  to  time  derivitivies  or  finite- 
differences. 


To  construct  the  projection  operator  using  product  of  linear  approximations  one 
can  define  two-dimensional  function  F(h,T)  as  the  following  sum  [25]: 

z” )  —  ^00  ^01^  ^10 

This  formula  defines  the  function  inside  an  arbitrary  (m,n)  rectangle  through  the 
values  of  the  function  at  the  four  angular  points  (x,y)=  {(0,0);  (0,1);  (1,0);  (1,1)}, 
where  x,  y  are  normalized  coordinates  x  =  (t  -  t^)  /  AT,y  =  (h-  h^)  /  Ah. 

The  next  example  of  projection  operator  design  is  based  on  the  product  of  local 
cubic  splines  [25].  Then  the  function  F(h,T)  takes  the  form 

F{h,T)=  (21) 


As  it  was  described  above  the  function  represented  through  the  normalized 
coordinates  x  =  (r  -  x^)  /  Ax,y  =  (h-  h„)  /  Ah,  inside  an  arbitrary  (m,n) 
rectangle,  can  be  obtained  through  the  values  of  the  function  at  the  four  angular 
points  (x,y)=  {(0,0);  (0,1);  (1,0);  (1,1)}: 

F(x,y)  =  +  F^P^i  +F^P^  +  FSPS  +  PA+. .... 

Here  F^ ,  F^  are  the  values  of  the  function  partial  derivatives  with  respect  to  x,y, 

F^  is  the  value  of  the  function  partial  derivative  of  the  second  order  with  respect 

to  X  and  y.  The  total  sum  (21)  will  contain  16  summands,  P^pi^^y)  are 

corresponding  polynomials  of  the  power  up  to  three.  Two  exaples  of  coefficient 
calculation  are 


Poo  =  4x -  ^x^y^  -  6x^y^  +  9x^y^  +  2y^  -  3y^  +  2x^  -  3x^  + 1 

Pqo  =  2x^y^  -  4x^y^  -  3x^y^  +  6x^y^  +  2xy^  -  +x^  -  2x^  +  x 

Then,  integrating  in  each  cell  of  the  given  polynomials  we  can  produce  the 
corresponding  elements  of  the  matrix,  as  in  the  case  (16-20).  The  cubic  spline 
approximation  is  differ  from  the  described  above.  Such  representation  makes  it 
possible  to  find  not  only  the  function  but  also  its  first  derivatives. 
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Now  we  have  describe  four  methods  of  the  RT  projection  operators 
constructing.  As  it  was  mentioned  above,  the  projection  operators  with 
approximations  of  higher  orders  allow  a  better  approximation  of  the  operator  of 
the  direct  probiem,  i.e.  they  enabie  us  to  come  closer  to  the  true  operator  of 
the  direct  problem  [25]. 

Description  of  the  program  for  caicuiation  of  the  different  versions  of  the  RT 
matrices  (operators)  for  the  high-orbit  RT  measurements 

2.  Program  <mat_rth1.for> 

This  program  designs  different  versions  of  the  operators  (matrices)  for  the 
TEC  and  TEC-difference  high-orbit  RT  measurements. 

Input  parameters  and  files: 

When  starting  program,  you  may  choose  the  type  of  RT  measurements 
(TEC  or  TEC-difference)  and  the  order  of  approximation  (one  of  four  mentioned  above 
approximation  methods). 

Fiie  <task.dir>  contains  the  same  parameters  as  the  input  file  for  the  dir_rth1.for 
program.  It  was  described  in  chapter  1 . 

File  <name_F.mat>  contains  names  of  output  files  with  parameters  of 
matrix,  nonzero  matrix  items  and  they  positions  in  matrix.  All  this  files  are  need  as 
input  files  for  program  of  the  inverse  problem  solution.  The  example  of 
<name_F.mat>  file  for  two-receiver  configuration  is  enclosed. 

Output  files: 

Files  <F_matr>  -  arrays  of  matric  of  corresponding  approximation 
for  each  receiver 

File  <Fparam>  -  parameters  of  matrices  for  each  receiver 
Files  <FJnt>  -  results  of  multiplications  of  calculated  matric 
and  model  structure 

File  <F_st>  -  array  (LST)  of  number  of  all  rays  which  cross  the 
corresponding  discret  of  modei  (SIRT  algorithm) 

Compilation 

mf77  mat_rth1.for  -ol  -486 
RUN 

ndprun  mat_rth1  .Itl 

1 .4.  The  solution  of  the  inverse  problem  for  satellite  RT  using  middie  and 
high-orbitai  satellites  of  opportunity. 


For  futher  computer  modeling  of  the  RT  problems  using  middle  and  high- 
orbital  satellites  of  opportunity  it  is  necessary  to  use  a  set  of  appropriate  electron 
density  distributions  models  of  the  protonosphere.  We  have  less  information  about 
structure  and  dinamic  of  magnetosphere  and  protonosphere.  Now  there  are  no 
detailed  and  precise  models  of  magnetosphere  and  protonosphere,  however  as  for 
ionosphere  ,  too.  But  there  is  no  need  in  it  for  estimation  of  possibilities  of 
radiotomography  of  the  magnetosphere  and  protonosphere.  To  solve  this  problem  the 
simple  models  of  magnetosphere  and  protonosphere  are  required,  which  include  the 
main  structural  features  (  plasmasphere  and  localized  natural  or  man-made 
irregularities  and  groups  of  irregularities  ).  However,  the  presented  package  of 
programs  makes  it  easy  to  design  many  other  structural  types  and  to  extend  this 
"ZOO"  as  far  as  possible  using  the  available  "detailes".  Here  we  will  use  the  simple 
protonospheric  model  ,  which  is  based  on  the  idea,  that  electron  density  distribution  is 
determined  mainly  by  geomagnetic  L  shells.  Fig.  14  illustrates  a  cutaway  view  of 
geomagnetic  L  shells  from  the  northen  mid -latitude  station  [23].  The  L  shells  greater 
than  4  are  considered  to  be  open  field  lines,  and  consequently  almost  do  not  contents 
the  electrons.  In  addition  supposed  is  the  possibility  of  existence  of  localized  natural 
or  man-made  irregularities  and  groups  of  irregularities.  In  coordinate  system  T  is 
counted  from  horisontal  ^is  -  counterclockwise),  such  electron  density  distribution, 
set  by  the  state  L=4  shells,  may  be  assigned  by  the  following  model  function: 

AT  =  iV„  +  iV,(l  -  /!  /  A,)  +  JV,  cos^(|  +  ^))  (22) 

First  two  items  describe  the  background  density  (interplanet  plasma,  solar  wind,  etc.) 
and  the  last  item  describes  the  plasmosphere.  In  computer  modelling  we  will  suppose 
that  (  if  there  are  no  another  values) 

Nq  =  Wcm  W^cnf^‘,N2  =  5*Wcrn^‘,}\  -  20000^m; 

h2  =  20000A:w;  Tj  =  5000km.  Let  us  underline  one  more  time  that  precise  functional 

dependencies,  describing  the  structure  and  location  of  L  shells,  are  unimportant  for 
RT  modelling.  The  qualitative  congruence  of  reconstructed  "pictures"  is  quite 
sufficient.  The  electron  density  distribution,  setting  by  model  22  is  cited  in  fig.  15, 
where  h  and  T  are  measuring  in  km. 

Let  us  turn  to  the  description  of  computer  modelling  results.  We  wil  consider 
three  main  experimental  geometries,  which  were  chosen  from  data  of  the  section  1.1: 
I  -  three  middle-orbital  and  one  high-orbital  satellites  (fig.?),  II  -  one  middle-orbital 


Earth  Rodii  (km) 

tCulowoy  view  of  the  earth  olong  the  stotion  meridion,  showing  the 
magnetic  L  sheds 'dnd  the  station  elevotion  onglesi  both  north  ond 
south  of  the  .stotion.  Dashed  line  represents  GPS  orbitol  heightens]. 
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Linear  RT,  Real  Georaetry  (h,tdu) 

Hodel  area:  K-f-SeBe ^,13000,0],  y-r0^,1?20a^J,  Grid  44x39  -  1716  itei»» 
Receivers  at  altitude  0^,  orbit  height  19200.0 
Mreceiver  Wrays  Positions  Orbit 

1  5  0  J  — >  i3223JB  1930.0  — >  3861.0 

2  5  386QJ3  — >  12674^  1938^  — >  3217  J 

3  5  7719  JB  — >  11685.0  1930.0  — >  2509.0 


and  two  high-orbital  satellites  (as  first  two  sateiiites  in  fig.6),  III  -  one  high-  orbital  and 
three  middle-orbital  satellites,  roteiting  in  opposite  direction  (fig.9).  Other  experimental 
schemes  are  either  unusefuii  for  reconstruction  of  two-dimensionai  crosses  ( 
contribution  of  high-orbitai  sateiiites  and  onearth  receivers),  or  give  the  same  results. 

RT  reconstruction  of  plasmaspheric  structure  or  location  of  L  shells  boundaries 
gives  quite  good  results  for  all  three  schemes  under  the  condition,  that  we  have  a 
priori  information  about  eiectron  density  magnitude,  location  and  sizes  of  desired 
structure.  However,  here  the  sufficient  accuracy  of  this  a  priori  information  is  required 
to  be  30-50%.  According  to  the  comparison  of  TEC,  restored  by  reconstruction  and 
real  TEC,  it  is  possible  to  chose  appropriate  initial  approximations  and  finite  results, 
varying  different  initial  approximations.  Fig.  16- 17  show  the  results  of  RT  reconstruction 
model  (22)  according  to  scheme  1,  where  the  Initial  approximation  contain  the 
mistakes  In  determination  of  real  distribution  parameters;  fig.  16  (the  mistake  in 
estimation  of  vertical  size  is  +10%),  fig.  17  (the  mistake  in  estimation  of  the  size 
is  -  15%),  fig.  18  (the  mistake  in  estimation  of  electron  density  maximum  is 
+15%).  The  experimental  scheme  was  made  in  following  way:  the  area  of  positive  T 
can  be  seen,  that  is  why  we  can  observe  the  only  initial  approximation  from  the  left  of 
reconstruction. 

The  quality  of  RT  reconstructions  of  smaller  structures  with  localized 
irregularities  is  notably  worse  according  to  scheme  I.  The  results  of  RT  reconstruction 
model  (22)  with  additional  localized  irregularities  are  cited  in  fig.  19-27.  The  localized 
irregularity  of  the  model  from  fig.  19  under  reconstruction  (fig.20)  is  extracting  and 
distoring.  The  irregularities  located  on  the  border  of  L=4  shell,  like  model  from  fig.21 , 
are  reconstructed  badly.  The  results  of  RT  plasmaspheric  reconstruction  with  such  a 
irregularity  with  different  magnitudes  of  local  electron  density  maximum  of 
irregularity  are  cited  in  fig.22  =  5*10’ ),  fig.23  =  10*10’ m“^),  fig.24 
{N^  =  25*10’ /w"^).  It  is  seen  that  even  the  irregularities  of  quite  high  concentration 
exceeding  maximum  inside  L=4  shell  in  2-5  times,  practically  can  not  be 
reconstructed.  Under  the  reconstruction  ,  the  local  irregularity  is  only  "spreading" 
along  the  enclosing  space.  The  same  picture  is  observed  under  the  reconstruction  of 
more  complex  model  with  three  local  irregularities  (fig.25,26).  Here  also  (fig.27)  the 
"Internal"  irregularities  practically  are  not  reconstructed,  and  the  "external"  are 
distorting  under  the  reconstruction.  Such  a  lower  quality  of  reconstructions  is 
connected  with  bad  experimental  geometry  (fig. 7),  which  is  cited  in  coordinates  (h,x) 
in  fig. 28.  The  crossing  angles  of  probing  rays  are  small  in  this  scheme,  what  leads  to 
strong  expansion  and  deformation  of  localized  objects  under  the  reconsruction. 


The  scheme  II  also  leads  to  the  "spreading"  of  reconstructed  objects.  It  can 
be  seen  well,  if  we  cite  an  example  of  restoration  of  model  structure  (22)  with 
localized  irregularity  in  its  center  (fig. 29),  which  practically  is  not  selected  being 
reconstructed.  The  same  pictire  of  "spreading"  can  be  observed  on  fig.31-32  under 
the  reconstruction  of  a  pair  of  localized  irregularities.  However  there  is  another  cause 
of  such  "spreading"  for  experimental  geometry  II,  than  it  was  in  I  geometry.  In  fact  it 
can  be  seen  from  experimental  scheme  in  coordinates  (/i,r),  that  the  cross  angles  of 
probing  rays  are  varying  in  a  wide  range,  but  ,  unfortunately,  there  (  in  this  scheme) 
exist  no  extensive  areas  of  uncrossed  rays,  which  contain  the  rays  of  the  only  one 
satellite.  The  existence  of  such  areas  leads,  while  reconstruction,  to  the  "spreading" 
of  restored  object  from  the  crossing  ray  area  to  the  area  of  uncrossed  rays.  It  is  easy 
to  understand,  that  if  there  exist  the  localized  irregularity  (shaded  in  fig.33)  in  the 
crossing  ray  area,  than  the  algorithms  will  localize  and  select  this  irregularity  correctly 
in  the  crossing  ray  area  under  the  reconstruction,  but  will  leave  its  traces  (  shaded 
horisontally  in  fig.33)  in  the  uncrossing  ray  area.  To  illustrate  this  thesis  let  us  consider 
the  simple  model  of  one  localized  irregularity  (fig.34),  the  result  of  its  reconstruction  is 
cited  in  fig.35,  where  extensive  traces  of  such  a  "spreading"  can  be  seen  from  the 
crossing  ray  area  to  uncrossing  ray  area.  However  such  a  "spreading"  can  be  avoided 
if  we  have  a  priori  information  about  negligibly  small  electron  density  in  the  uncrossing 
ray  area,  that  is  the  irregularities  are  located  only  in  the  crossing  probing  ray  area. 
The  example  of  RT  reconstruction  of  tha  same  model  (fig.34)  with  the  use  of  such  a 
priori  information  is  cited  in  fig.36. 

The  results  of  RT  reconstructions  for  scheme  III  with  satellites,  rotating  in  the 
opposite  directions,  and  with  quasifan  beams  are  little  bit  better.  The  odd  figures 
(fig.37-48)  show  the  protonospheric  models  with  various  localized  irregularities,  even 
figures  show  the  results  of  RT  reconstructions.  Here,  as  usual,  the  localized 
irregularities  are  selected  under  the  reconstruction  though  the  distortions  and 
attefacts  are  significant.  The  higher  quality  of  RT  reconstructions  of  scheme  III  is 
explained  by  the  better  experimental  geometry  (fig.9,  cited  in  fig.33,  in  coordinates 
(h,r),  from  which  the  large  crossing  angles  of  probing  rays  and  small  uncrossing  ray 
areas  on  the  border  of  reconstruction  area  can  be  seen.  Surely,  large  ray  curvature 
on  borders  of  reconstruction  area,  the  existence  of  upper  and  lower  ray  beams  lead 
to  corresponded  distortions  and  artefacts  in  RT  reconstructions.  However  the 
development  of  effective  methods  and  algorithms  of  RT  reconstructions  in  connection 
with  the  similar  certain  schemes  makes  it  possible  to  improve  the  results. 


Linear  RT,  Real  Geonetry  (h,tau) 

Model  area:  1-9000 .0,9300^],  9-10^,18100^1,  Grid  44x39  -  1716  ibe»s 
Receivers  at  altitude  6.0,  orbit  height  18100.0 
Hreceiver  Hrays  Positions  Orbit 

1  7  ^8815.0 —>  6611.0  -4266.0 —> -1865.0 

2  7  -7374.0— >8815.0  1000.0  —>  3519.0 
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Lineup  RT.  Resil  GeoMetry 

Model  area:  X-[-5508£,11^£],  V^lBJUJZtmah  6rid  44x39  -  1716  iteas 
Receivers  at  altitude  BJB,  orbit  height  21480^0 
Nreceiver  Hrays  Positions  Orbit 

1  7  365B£— '>-5814^  2434X1 —>  3537  JB 

2  7  7300^ —> -4637^  2434^ —>  3923^ 

3  7  109S0A  — > -2647A  2434^— >4889^ 


Description  of  the  program  for  soiution  inverse  problem  for  TEC  RT  or 

TEC-difference  for  high  orbital  RT  . 


Program  <slv_rth1.for> 

Input  parameters  and  files. 

When  starting  program,  you  may  choose  the  type  of  RT  measurements. 

Answer  =  1  :  TEC  RT 

Answer  =  2  :  TEC-difference  RT 

NF,  NR  -  number  of  discrets  on  the  horizontai  and  verticai  grid 
NREC  -  number  of  receivers 

’F_MOD'  -  input  fiie  of  model  structure  (program  <dir_rth1.for>) 

'F_XO'  -  input  fiie  of  initial  guess 

(program  <dir_rth1.for>  or  function  <guess>) 

'Fparam'  -  input  file  with  parameters  of  matrix 
XIST  -  model  structure  from  input  file  'F_MOD' 

File  'name_F.slv'  containes  Nrec  lines  with  names  of  input  and  output  fiies. 

The  example  of  this  fiie  is  enclosed. 

line  1  -  name  of  file  with  parameters  of  matrix  (Fparam) 

line  2:  namel  -  name  of  input  fiie  with  matrix 

name2  -  name  of  input  file  with  either  TEC  or  TEC-difference 
('Finput') 

names  -  name  of  output  file  with  results  of  multiplications 
of  input  matric  and  reconstructed  structure  ('Fout') 
line  Nrec-i-2:  name  of  file  with  number  of  aii  rays  which  cross 
the  corresponding  discret  of  model 
line  Nrec+3:  name  file  with  model,  name  file  with  guess, 
name  file  with  reconstruction 
line  Nrec+4:  name  file  with  errors 
Npi  -  array  of  contants  (2pn)  for  each  receiver 
er_n  -  levei  of  noise  in  % 

Npi,  er_n  are  input  from  the  screen 

Nsoive  -  the  method  of  solution,  it  is  input  from  the  screen 

REL  -  array  of  relaxation  parameter 

LMAX  -  max  number  of  nonzero  matric  elements  for  each  receiver 
AZ  -  array  of  nonzero  matric  elements  for  all  receivers 
NST  -  array  of  corresponding  column  nonzero  matric  eiements  for 
all  receivers 

LST  -  array  of  number  of  all  rays  which  cross  the 

corresponding  discret  of  model  (SIRT  algorithm)  from 
input  fiie  <F_ST>,  it  is  similar  to  LST  in  program  <mat_rth1  .for> 

Niter  -  number  of  iterations,  it  is  input  from  screen 
Subroutine  <DEFSYS>  -  solution  of  linear  system  equations  for 
one  ray. 

Subroutine  <ercl2>  calculates  errors  in  metric  C  and  L2 
Subroutine  <VMINM>  calculates  min  and  max  of  array 
Function  <RAN>  calculates  random  values  in  [0,1] 

Function  <GUESS>  -  for  initial  guess 

RMAX,  RM,  Zmax,  ZSM,  B1,  B2  -  parameters  for  function  <GUESS> 
Output  files: 

File  <F_REC>  -  file  with  the  reconstruction 
Fiie  <er_solv>  -  errors  of  right  items  and  reconstruction  in 
metric  C  and  L2 

Fiies  <Fout>  -  results  of  multiplications  of  the  matric 
and  reconstructed  structure 


Execution 

mf77  siv  rthl.for  -ol  -486 
RUN 

ndprun  slv_rth1.ltl 
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Part  2 

Algorithms  and  Programs  for  Nonlinear  Radiotomography. 


2.1.  The  Soiution  of  the  Direct  Problem  of  Radio  Wave  Propagation  for 

Nonlinear  Radiotomography. 

The  presence  of  strong  gradients  of  electron  concentration  leads  to  the  necessity 
of  using  the  ray  refraction  in  the  problems  of  ray  radiotomography.  The  ignoration  of 
refraction  can  give  wrong  results  of  RT  reconstructions,  and  also  limits  the  resolving 
power  of  RT  system.  Taking  the  refraction  into  account  it  transforms  to  more  complex 
problems  of  RT.  Supposed  that  the  set  of  integrals  over  retractable  rays,  which 
curvature  is  determined  by  propagating  medium  are  known. 

TEC=  \Nd(T=L[N]  (1) 

L[N] 


Where,  as  it  was  mentioned  earlier,  J  Nda  denotes  the  integration  of  the  electron 

density  N  along  the  signal  path.  However  in  general  case  the  ray  paths  L[N]  now 
depend  on  the  distribution  of  electron  concentration  and  are  determined  by  ray 
equations.  The  theory  and  methods  of  linear  RT  reconstructions  has  been  given  in 
our  earlier  papers  [1 -3,6,8].  Similar  tomographic  problems  on  linear  integrals  have 
been  solved  in  other  fields  of  science  and  the  methods  of  solving  such  problems  are 
well-known.  In  the  case  of  ionospheric  diagnostics,  however,  the  major  difficulty  is 
that  the  linear  integral  of  the  electron  concentration  is  proportional  to  the 
absolute  phase  O,  whose  accurate  determination  is  practically  impossible, 
therefore  we  have  suggested  [1-3,6,8,25]  to  apply  the  phase-difference  (or 
Doppler)  RT  approach,  which  has  given  quite  satisfactory  results.  Note  that  to  date 
most  works  on  the  ray  RT  have  been  based  on  the  phase  approach,  which  can  not 
provide  high-quality  and  reliable  results.  It  was  not  until  quite  recently  that  an 
acceptable  alternative  has  been  suggested,  that  of  reconstruction  by  the  relative 
phase  (relative  TEC)  data  [26,27].  Such  an  approach  Is  possible,  but  we  belive, 
however,  that  the  sensitivity  of  this  method  would  be  lower  than  that  of  the 
phase-difference  (or  Doppler)  RT  [25].  Such  an  approach  also  can  be  usefull  in  a 
case  of  nonlinear  RT.  Than  we  will  use  similar  system  of  equations  where  data  will  not 
be  TEC  or  absolute  phases  ,  but  difference  of  phases  or  Doppler  frequences  Dia^), 


which  are  proportional  to  derivative  of  TEC  with  the  respect  to  tiie  angle  of  satellite 
that  changes  lineary  for  cyclic  orbite  with  the  respect  to  time. 

In  this  article  we  will  consider  the  direct  problem  of  ray  path  calculation 
according  to  electron  concentration  N.  As  it  was  mentioned  earlier  we  will  use  the 
coordinate  system  (h,T).  To  calculate  ray  paths,  differential  equation  system  should 
be  used  [28-29].  Here  it  would  be  convenient  for  us  to  use  the  following  equations  for 
calculations: 

1.  Equations  that  calculate  the  ray  path  directly  with  the  respect  to  "angle  of 
refraction”^ . 


dr 


—tanO 

h 

i  f  ^  h  d{nh)  dr 
nh\  ^  R  ch  dhj 


(3) 


2.  "dinamic  equations"  ,  describing  the  ray  path  by  the  "time"  parameter  t,  which 
physical  meaning  is  the  optical  lenght. 


dh  {h  +  R)*T^ 

dt  ~  R}  ^  2  Sh 
dUh  +  Rf  '  dt 


T 


dt 


(4) 


Here  R  -  radius  of  the  Earth,  h  -  the  height  over  the  Earth's  surtece,  n  -  coefficient  of 


a 


refraction,  «  =  (if  N  is  given  in  units  10^^ m  f  -  in  MHz, 


a  »  80.611),  T  -  distance  along  the  Earth's  surface,  (remember  that  h,  r  are 
analog  of  polar  coordinates,  connected  with  the  center  of  the  Earth). 

To  estimate  the  area  of  using  of  this  methods,  which  help  to  draw  the  ray  path, 
test  and  comparison  of  both  methods  were  made.  To  model  the  ray  path,  the 
quasireal  models  of  ionosphere  N  (shown  in  fig.1)  were  used.  In  fig.1  shown  are  the 

through  and  local  irregularities  ,h  g  [0,1000],  r  g  [-1000,2200]  km.  Here  and  later  all 
the  pictures  in  shades  of  grey  are  illustrated  in  Cartesian  coordinates,  so  that 
coordinates  (h,T)  look  as  curvalinear  coordinates,  and  unrefractable  rays  are  linear. 


The  notes  of  corresponded  distances  over  h  and  x  in  the  pictures,  illustrating  the 
paths  ,  are  not  shown,  as  the  only  essential  things  here  are  the  distortion  and 
divergence  of  rays.  Let  us  cite  the  results  of  investigations  for  three  essentially 

different  cases  of  RT  probing  of  ionosphere  (concentration  is  cited  in  units  W^rn^, 
frequency  of  probing  f  in  MHz). 

1.  "Lower"  electron  concentration  (with  lower  and  average  solar  activity,  max 

N  =  0.5*10‘^m  =  mMHz). 

2.  "High"  electron  concentration  (with  high  solar  activity,  max 
N  =  5*W^nf\f  =  mMHz) 

3.  High  electron  density  max  N  =  and  lower  frequency  of  probing 

f=50  MHz. 

Charts  of  A/i(r), Ar(/i), A7!E'C(i:o), Aji!)(ro)  dependences  were  dedused  to 
illustrate  path  calculations.  Here,  A  -  deflection  of  reflected  ray  from  "direct"  ray  ( 
may  be  negative,  because  the  magnitude  for  "direct"  ray  was  subtracted  from  the 
corresponded  magnitude  for  real  case),  constructed  with  no  respect  to  nonlinear 
effects,  D-Doppler  data.  There  are  pictures  which  have  ray  paths  with  the  illustration 
of  ionosphere  in  grey  shades,  where  the  influence  of  inhomogeneities  and  parameters 
of  the  problem  (N,  f)  on  the  ray  path  is  seen. 

To  test  methods  of  ray  calculations  it  is  expedient  to  use  so  to  say  "quasiexact 
solution"  for  layered  ionosphere.  Easy  to  see,  that  in  a  case  of  spherically  layered 
ionosphere,  equations  (4)  can  be  integrated  and  it  is  possible  to  get  expressions  for 
ray  path  x(]i)  and  integral  I  over  h. 


Ti  -  T(h)  = 


cos  pdr 


r  -  cos^  P  ~  N(R))r 

V  J 


(5) 


R+hg 

1=  J 

R 


I 


-^[N(r)-N(R)]r^ 


r 


cc 

-—2  [N(r)  -  N(R)]r^  -  R^  cos^  p 


N(r)dr 


r 


The  calculation  of  such  "quasiexact  solution"  were  carried  out  for  quasireal  layered 
ionosphere  ( the  model  of  fig.1,  without  the  through  and  irregularities).  In  table  1  the 


values  of  maximum  deflection  of  rih^)  in  km  on  the  height  of  satellite  are  shown.  It 
was  found  by  means  of  both  methods  from  differential  equation  system  with 
"quasiexact  solution"  for  one  case  of  layered  ionosphere  and  at  various  angles  ^ 
from  vertical  line.  Moreover  the  deflection  over  r  has  the  multiplicity  of  100  m  at  fall 
angle  ^  =  80°  ,  what  is  quite  satisfactory  because  at  such  angies  the  ray  passes  the 
distance  multiplicity  of  3000km  in  ionosphere.  A  step  of  solution  of  differential 
equation  system  by  method  1  or  2  here  is  equal  to  100m  along  the  ray  lenght,  so  the 
deflection  has  the  multiplicity  of  1  step.  It  is  seen  that  the  second  method  gives  better 
results. 

To  compare  the  deflections  of  t  of  solutions  T(h),  which  were  received  by  two 
methods  "in  emptiness"  (n=1),  from  exact  solution  "in  emptiness"  (fig.2),  the  first 
method  ,  iliustrated  by  dotted  iine  in  picture,  led  to  maximum  deflection,  especially  at 
large  fall  angles  of  ray.  It  is  connected  with  the  fact,  that  the  step  by  step  solution  of 
differential  equation  system  was  accompanied  many  times  by  operation  tan.  This 
leads  to  essential  supply  of  mistakes  because  of  the  finite  bit-word  lenght  of  computer 
machine.  And  though  the  second  method  doesn't  aliow  to  chose  fixed  step  over  h  and 
,  with  the  solution  of  differential  equation  system,  it  doesn't  lead  to  multiple 
calculations  of  trigonometrical  functions,  and  that  is  why  it  is  more  suitable  for  the 
solution  of  RT  nonlinear  problems  from  the  point  of  view  of  accuracy  and  speed. 

As  it  was  expected  in  the  case  of  lower  electron  concentration  at  maximum 
0.5*10*^ and  f=150  MHz  (case  1)  the  influence  of  ionosphere  on  ray  path  is 
insignificant  (fig.1).  Fig. 3-6  show  the  charts  of  the  dependences  Ah(T),AT(h)  for  fall 
angles  ^  =  40,80°  from  vertical  line  (pic.3,5  Ahit)  for  ^  =  40,80°  correspondly, 
pic.4,6  Arih)  for  ^  =  40,80°  correspondly).  Pic.7  shows  the  values  of 
TEC{tq),TECq{tq),  and  fig.8  shows  D(tq),Dq(tq),  where  TECq(tq),Dq(tq)  - 
magnitudes  of  TEC  and  Doppler  for  "direct"  ray  (illustrated  by  dotted  line).  It  is  seen 
that  with  the  increasing  of  the  fall  angle  (  what  is  corresponded  to  greater  distance  of 
ray  in  ionosphere),  deflection  of  ray  from  the  "direct"  one  over  t  increases  (compare 
fig.4  and  6),  and,  correspondly,  the  difference  between  magnitudes  of  TEC  increases. 
But  even  for  ^  =  80°  the  difference  between  TEC(tq)  and  TECq(tq)  doesn't 
exceed  0,1%  for  this  case,  and  At  doesn't  achieve  2km  (fig.6). 

In  the  case  of  high  electron  concentration  {N  =  =  ISOMHz  -  the 

second  case)  the  influence  of  nonlinear  effects  on  ray  path  is  essentially  stronger.  Fig. 
9-12  show  the  charts  of  dependences  Ah(T),Ar(h)  for  fall  angles  ^  =  40,80°  from 
vertical  line  (fig.9,11-  for  Ah(T)  for  ^  =  40,80°  correspondly,  fig. 10,12  At(/i)  for 
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)opler,  rad/sec  TEC, 
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^  =  40,80°  correspondly).  Fig. 14  shows  the  values  of  TEC(tq)  and  TECq(tq),  and 
pic.  15  -  the  same  for  Z)(z-o),2)o(ro).  In  addition  the  mistake  in  data  (TEC)  at  a  fall 
angle  equal  to  80  degrees  is  now  essential  and  it  reaches  3%  (fig.  13)  and  Ar(A)  is 
attaining  to  18  km  (fig.  14). 

The  estimations  of  small  deflections  from  ray  path  in  emptiness  can  be 
received  without  any  difficulties  from  equations  (4).  Let  us  consider  the  vertical  ray  fall. 
It  can  be  shown  thgt  the  estimations  of  angles  of  vertical  deflections  and  deflections 
themselves  would  be  the  same  for  another  direction  of  ray  fall.  The  angle  of 

vertical  deflection  for  normal  ray  fall  x  —  4Pz-/  P  expressed  by  the  analog  of 

impulse  for  "dinamic"  equations  (4). 

Let  us  remember  that  =  h,  =  (f  /  ^ —  Ph  +  Pt  ^  P 

and  for  quasivertical  ray/^  fn  Pfj(t  =  0)  =  1 .  The  increase  of  impulse  according  to 

T  is  determined  by  second  equation  (with  the  respect  to  ignoration  of  difference 
between  (R+h)/R  and  1 ). 


Az^AP^ 


1 

2  dr 


At 


2  dr 


a 


a  AN 


aAN 


a 


a 


a 


2f^  dx  2f^  a  /  2 


(6) 


Supposed  that  At  oi  a  (because  dt  =  dl  /  n  ^  M),  the  size  of 

inhomogeneity,  and  concentration  gradient  is  equal  to  ratio  of  its  overfall  AN  and  the 
half  of  total  size  of  inhomogemeity.  Then  the  transverse  ray  deflection  (distance 
between  rays)  would  be  propotional  to  the  passed  distance  L. 


Ax  ^ 


-^ANL 

P 


(7) 


What  is  (for  probing  frequency  150  MHz)  equal  to: 


Ax  «  3,6*10~^ANL 


(8) 
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Here  AN  is  measuring  in  units  of  10  m  .  It  is  notable  that  the  relations 
(6,7)  do  not  involve  the  sizes  of  irregularity,  everything  is  determined  by  overfall  of 

12 

electron  concentation.  For  typical  overfalls  AN  ^  (1  —  3)10  and  ray  lenght 


L  (1  —  3)10  km  the  deflection  (8)  would  be  varyed  In  the  range  between  one 
and  tens  of  kilometers. 

Fig.  15  shows  the  deflection  AT(h)  (when  ^  =  80°),  calculated  for  the  same 
layer  (fig.1)  in  the  third  case  of  high  concentration  and  lower  frequency  probing.  For 

the  second  and  the  third  cases  the  parameter  aMl  /  f  differs  in  9  times,  that  is 
why  the  deflection  ArCh)  behaves  like  the  deflection  from  fig.6,  but  nearly  ten 
times  more.  In  the  scale  fig.1  the  deflections  equal  to  10  km  ,  which  are  typical  for  the 
second  case,  are  unnoticable,  that  is  why  to  illustrate  path  calculations  the  figures  for 
the  third  case  are  cited.  Fig.  16  shows  well  the  "focusing"  of  ray  by  the  electron 
concentration  trough  .  Rg.17-18  show,  that  the  local  irregularities  can  deflect  the  ray 
in  various  directions  from  the  straight  ray,  and  it  depends  on  from  what  side  relatively 
to  irregularity  maximum  the  ray  passes. 

As  it  follows  from  the  results  of  this  section,  the  small  angular  and  linear 
deflections  of  real  rays  from  linear  approximation  are  propotional  to  variations  A^of 
electron  concentration  and  do  not  depend  on  the  size  of  irregularities.  For  probing 
frequancy  150  MHz,  these  deflections  are  estimated  less  than  a  kilometer  for 
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variations  AN  <  10  m  ,  that  is  when  the  solar  activity  is  not  high.  For  example, 

12 

these  deflections  can  reach  tens  of  km  for  AN  >  (2  —  3)10  m  ,  when  the  solar 
activity  is  high.  Hence,  it  follows  that  for  RT  investigations  on  150  MHz  the  respect  of 
refractable  corrections  is  necessary  only  for  high  solar  activity  and  for  the  sizes  of 
digitization  elements  less  than  20-30  km.  The  use  of  lower  frequency  radiowaves 
f  <  (50  —  100) MHz  for  RT  investigations  requires  the  respect  of  refraction. 


Table  1. 
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METHOD  1 
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METHOD  2 

0 

3.78E-5 

3.78E-5 

20 

3.08E-2 

4.47E-3 

40 

-7.14E-2 

1.17E-2 

60 

-0.12 

3.03E-2 

80 

-0.27 
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Description  of  the  programs  for  calculation  TEC,  Phase  and  Phase-difference 

measurements  for  Nonlinear  RT. 


System  Requirements 


-  Computer:  IBM  AT-486  or  compatible  (with  coprocessor) 

-  Operating  System:  MS-DOS  or  PC-DOS  version  3.0  and  later 

-  Memory:  at  least  Extended  memory  8  Mbytes 

(depends  on  geometry  and  type  of 
approximation  reconstructed  function) 

-  Hard  Disk  Space:  35  Mbytes 

-  Software:  NDP-FORTRAN-486  Compiler  Ver.3.1 

NDPUnk  Ver.3.0  (C)  MicroWay,  Inc. 


Program  <dir_nl.for> 

This  program  solves  direct  problem  for  Nonlinear  and  Linear  RT  measurements, 
namely,  determines  the  model  structure  and  calculates  the  TEC,  Phase  and 
Phase-difference  integrals  on  the  model  structure. 

Input  parameters  and  files: 

When  starting  program,  you  may  choose  the  type  of  RT  (Nonlinear  or 
Linear)  and  RT  measurements  (TEC,  Phase  or  Phase-difference).  This 
program  calculates  Phase-difference  measurements  only  for  linear  RT.  Then 
program  asks  about  integrals  to  calculate.  If  you  need  only  model  structure  you  can 
answer  'No’. 

Input  file  <task.dir>  contains  information  about  RT  measurements  geometry 
and  defines  the  names  for  output  files.  The  example  of  this  file  is  enclosed.  It 
contains  the  following  lines: 

'F_Nmod'  -  input  file  that  describes  parameters  of  model  structure. 

The  example  of  this  file  is  enclosed. 

’F_moder  -  output  file  with  model  structure  (GRD-format) 

NF,  NR  -  number  of  discrete  on  the  horizontal  and  vertical  grid 
Nrec  -  number  of  receivers 

FP,  TSATO,  TSAT1,  RAY  -  Nrec  lines  with  coordinates  of  receivers, 
initial  and  final  positions  of  satellite 
in  km  and  number  of  rays  for  each  receiver 
HFIST,  HFIN  -  vertical  size  of  reconstructed  area  in  km 
TFIST,  TFIN  -  horizontal  size  of  reconstructed  area  in  km 
NJ  -  number  of  the  discrete  steps  along  the  ray 


TINT’  -  Nrec  lines  with  the  names  of  output  files  for  integrals 
for  each  receiver 

Freq  -  the  probing  frequency  (MHz) 

ts,  ddt,  ddh  -  parameters  for  subroutine  <INTERG> 

The  model  structure  is  determined  by  function  <FMODEL>  which  uses 
functions  <FUNC1>,  <FUNC2>,  <HOMPAR>,  <HOMCOS>.  The  combinations 
of  these  functions  make  it  possible  to  obtain  different  model 
structures. 

Subroutine  <DEFPSI>  determinates  the  angles  of  rays  on  the 
satellite's  cooordinates. 

Subroutine  <DEFINT>  calculates  the  integrals  for  all  rays  for  each 
receiver  and  uses  subroutine  <INTERG>  which  calculates  the 
integral  for  one  ray. 

Output  files: 

MODEL.GRD  -  file  with  model  structure  (GRD-format) 

RIes  <linjnt>  -  arrays  of  either  TEC,  Phase  or  Phase-difference 
integrals 

Compilation 

mf77  dir_ni.for  -ol  -486 

RUN 

ndprun  dir_nl.ltl 
Program  <dpl_ni.for> 

This  program  calculates  Phase-difference  measurements  for 
nonlinear  RT  for  one  receiver. 

Input  parameters  and  files. 

'SFILE'  -  input  file  with  phase  array 
NZ  -  dimension  of  phase  array 
Subroutine  <rd>  reads  of  input  file 

Subroutine  <dpl>  calculates  doppler  by  means  of  approximation  of 
phase  of  polynomial  of  the  power  3  and  uses  subroutines  <regres> 
and  <gauss>. 

Subroutine  <wr>  writes  in  output  file. 


Output  file. 


'NFILE'  -  output  file  with  doppler 


Compilation 

mf77  dpl_nl.for  -ol  -486 
RUN 

ndprun  dpLnI.ltl 
Program  <tblint.for> 

This  program  calculates  linear  interpolation  of  input  array, 
input  file. 

’interp.ts'  -  input  file  contains  number  of  receivers  and 

Nrec  lines  with  names  of  input  and  output  files, 
with  dimensions  of  input  and  output  arrays. 

The  example  of  this  file  is  enclosed. 

Output  file. 

'F_appr'  -  output  files  with  interpolated  arrays. 

Compilation 

mf77  tbiint.for  -ol  -486 

RUN 

ndprun  tbiint.ltl 


2.2.  Design  of  the  different  versions  of  the  RT  operators  for  noniinear 

Radiotomography. 

The  problems  of  the  design  of  the  different  versions  of  the  RJ  operators 
(matrices)  are  considered  in  this  section.  At  first  it  is  necessary  to  perform 
digitization  of  the  linear  integrals  (1.1)  according  to  the  position  of  the  satellite, 
which  is  dependent  on  the  coordinates  or  the  angle  oto;  =  ■  The 

dependences  r(/i)  or  hij)  are  determined  from  ray  equations  (4).  A  series  of 
elevation  of  the  satellite  for  i-receiver  are  determined  from  relation  (1.6)  which 
define  a  series  of  discrete  values  of  the  linear  integrals  /,y  = The 


discretization  procedure  is  fully  described  in  [3,  25],  As  it  was  shown  earlier  in 
[3,. 8,. 25],  analogously  we  obtain  the  system  of  linear  equations,  which  may  be  either 
overdetermined  or  sub-definite. 


_  T 

^ij  ^  m/i  ^  ij  9 


or  ^MJ^M  -  Ij' 


(9) 


Here,  "renumbering"  of  the  ray  (i,j)->J  and  the  sells  of  the  ionosphere  (m,n)->M 
is  performed  in  the  second  equation.  The  problem  of  ionospheric  RT  according  to 
phase-difference  or  Doppler  measurements  require  higher-order  interpolation  than 
the  piecewise-constant  representation  of  the  function.  The  Doppler  frequency 
/  dt  measured  in  the  experiment  is  determined  by  the  phase  derivative  (2), 

namely: 


o  -  ^0 

i?  +  /jo  dt 

where  F,, -velocity  of  a  satellite  moving  uniformly  along  a  circular  orbit.  As  it  was  shown 
in  [3,.  10],  if  the  matrix  of  the  direct  problem  ->  I j  is  continuous  with 

respect  to  the  angle  of  the  satellite  a^,  it  is  possible  to  obtain  a  system  for 
phase-difference  or  Doppler  data  by  differentiating  (9)  with  respect  to  the  angle 
ao; 


Here,  Dj  =  AIj  /  Aa^  are  Doppler  data  and  Aj^f  =  ALj^^  /  Aa^  is  finite  -difference 
ratio  (or  derivative)  of  the  matrix  Lj^f  to  the  increment  of  the  angle. 

The  methods  of  constructing  the  projection  operators  or  Ljj^  matrices  for  the 
phase  RT  and  also  for  the  phase-difference  RT  are  similar  to  the  methods  presented 
in  details  earlier  in  [3,  25],  but  it  is  necessary  to  take  into  account  the 
dependence  T{h)  or  A(r).  The  following  operators  (matrices)  for  solving  the  RT 
problems  are  described: 

A  -  is  built  with  the  piece-constant  approximation, 

B  -  is  built  with  tiie  piece-planar  approximation, 

C  -  is  built  with  the  linear  product  approximation. 


D  -  is  built  with  the  cubic  spline  product  approximation. 

Projection  operators  with  approximations  of  higher  orders  allow  a  better 
approximation  of  the  operator  of  the  direct  problem,  i.e.  they  enable  us  to  come 
closer  to  the  true  operator  of  the  direct  problem.  It  is  necessary  to  note  that 
formally  in  the  shape  of  equations  (9-10),  the  problems  of  linear  and  nonlinear 
RT  are  identically.  But  it  is  not  so,  when  designing  of  operators  in  the  problems  (9- 
10)  we  use  direct  rays  in  the  linear  RT,  on  the  other  hand  in  the  nonlinear  RT  we  use 
the  rays  are  determined  by  system  (4).  Therefore  the  equations  (9-10)  relatively 
sought  function  F  (the  distribution  of  the  electron  density)  are  nonlinear,  because 
the  distribution  of  the  electron  density  defines  the  structure  of  operators  A  and  L  in 
(9-10).  Of  course,  if  the  curvature  of  the  rays  is  not  great,  than  the  operators  of 
linear  and  nonlinear  RT  will  be  practicallly  coincided. 

Table  1,2  show  examples  of  errors  of  calculations  of  the  phase  (integrals- 
Ij  )  for  each  of  receiving  site  (in  this  case  a  chain  of  four  receiving  sites  with 
cooordinates  r:  0,  300,  800  and  1200  km  are  used)  with  refraction  and  without 
refraction  for  model  1 :  two  irregularities  are  on  the  edge  of  the  trough  (this  model 
will  be  discribed  in  the  next  section)  with  the  help  of  different  operators:  B,  C.  For 
model  1  the  electron  density  AN  changes  from  1*10^^  (tabl.1)  to  5*10^^ 
(tabi.2)  at  probing  frequency  150  MHz.  Errors  of  the  numerical  simulation  can 
appropriately  be  characterized  by  number  d ,  which  shows  the  deviation  of  the 

calculated  function  (being  reconstructed)  F  from  the  true  function  F: 
5  =  ||f-?||/||f1i.  The  norms  of  the  spaces  F  and  /”  can  be  helpfully  used 

-  S(l^)  and  =  5(r)). 

As  it  is  seen  from  tabi.1  if  AN  is  about  10^  the  errors  of  caicuiation  of  the 
direct  problem  operator  with  and  without  refraction  are  commensurately  and  consist 
of  part  of  percent.  In  other  words  true  phase  little  differs  from  the  phase  calculated 
without  refraction,  that  is  illustrated  in  Fig.  19  (first  receiver)  and  Fig. 20  (third 
receiver).  The  true  phase  is  shown  by  solid  line,  and  dashed  line  is  the  phase 
calculated  without  refraction  in  these  figures.  The  curves  are  practically  coincided. 
The  phase  calculated  with  refraction  is  not  presented,  because  it  is 
indistinguishably  from  true  phase. 

We  have  another  situation  if  AN  is  about  5*10^^  As  it  is  shown  in  tabl.2,  the 
errors  of  calculation  of  the  direct  problem  operator  without  refraction  are  25%  in 

space  r  at  the  same  time  the  refraction  corrections  allow  to  decrease  these  errors 
up  to  units  of  percent.  The  curves  of  true  phase  (solid  line),  phase  calculated  with 
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refraction  (dotted  line)  and  phase  calculated  without  refraction  (dashed  line)  are 
presented  in  Fig.21 -24  for  all  four  receiving  sites.  One  can  see  essential  difference 
between  true  phase  and  phase  calculated  without  refraction.  In  particular,  the  phase 
calculated  without  refraction  has  different  behaviour  than  true  phase  in  Fig.21 
(receiver  1)  and  in  Fig. 24  (receiver  4),  and  the  difference  of  values  is  about  700-800 
radian  (Rg.  24).  In  other  hand,  the  phase  calculated  with  refraction  has  not  essential 
differences  from  true  phase.  Hence,  as  it  was  mentioned  in  section  2.1,  it  is 
necessary  to  take  into  account  the  refraction  corrections  for  RT  at  probing 
frequency  150  MHz  >A4ien  AN  >2-  i.e.  in  years  of  high  solar  activity. 


Table  1 

Accuracy  of  construction  of  the  direct  probiem  operator 


Model  1,  AN  « 


Nonlinear  RT  Linear  RT 


Operator 

type 

Receiver 

^2 

^2 

B 

1 

0.0084 

0.0061 

0.0131 

0.0096 

2 

0.0079 

0.0057 

0.0111 

0.0081 

3 

0.0068 

0.0051 

0.0123 

0.0073 

4 

0.0070 

0.0052 

0.0134 

0.0087 

C 

1 

0.0071 

0.0057 

0.0130 

0.0093 

2 

0.0081 

0.0059 

0.0110 

0.0082 

3 

0.0067 

0.0051 

0.0121 

0.0071 

4 

0.0069 

0.0051 

0.0115 

0.0086 

Tabie  2 

Accuracy  of  construction  of  the  direct  problem  operator 


Model  1,  AN  « 10‘^m  ^ 


Nonlinear  RT  Linear  RT 


Operator 

type 

Receiver 

^2 

^2 

B 

1 

0.0115 

0.0081 

0.2056 

0.1192 

2 

0.0098 

0.0070 

0.1486 

0.0769 

3 

0.0106 

0.0061 

0.1241 

0.0683 

4 

0.0135 

0.0065 

0.2591 

0.1242 

C 

1 

0.0092 

0.0073 

0.2041 

0.1191 

2 

0.0101 

0.0072 

0.1471 

0.0767 

3 

0.0099 

0.0060 

0.1219 

0.0679 

4 

0.0109 

0.0064 

0.2558 

0.1238 

Description  of  the  program  for  calculation  of  the  different  versions  of 
the  matrix  Coperators)  for  Nonlinear  RT 

System  Requirements 


-  Computer:  IBM  AT-486  or  compatible  (with  coprocessor) 

-  Operating  System:  MS-DOS  or  PC-DOS  version  3.0  and  later 

-  Memory:  at  least  Extended  memory  8  Mbytes 

(depends  on  geometry  and  type  of 
approximation  reconstructed  function) 

-  Hard  Disk  Space:  35  Mbytes 

-  Software:  NDP-FORTRAN-486  Compiler  Ver.3. 1 

NDPLink  Ver.3.0  (C)  MicroWay,  Inc. 

Program  <mat_nl.for> 

This  program  designs  different  versions  of  the  operators  (matrics)  for  phase 
RT  or  phase-difference  RT. 


Input  parameters  and  files: 

When  starting  program,  you  may  choose  the  type  of  RT  (Nonlinear 
Linear)  and  RT  measurements  (TEC,  Phase  and  Phase-difference). 

Input  file  <task.dir>  is  same  as  in  program  <DIR_NL.FOR> 

Parameters  <NF,  NR,  NREC,  NJ,  FP,  TSATO,  TSAT1,  RAY,  HFIST,  HFIN, 

TFIST,  TFIN,  Freq,  ts,  ddt,  ddh>  are  similar  to  same  parameters  of 
program  <DIR_NL.FOR>. 

nie  <name_F.mat>  contains  NREC  lines  with  names  of  output  files. 

'F_moder  -  input  file  of  model  structure  (program  <DIR_NL.FOR>) 

APTYPE  -  type  of  approximation  of  reconstructed  structure,  it  is 
introduced  from  the  screen 

Subroutine  <DEFPSI>  determinates  the  angles  of  rays  on  the 
satellite's  coordinates. 

Subroutine  <DEFMAT>  determinates  the  matrix  elements  with 
corresponding  approximation  for  one  receiver  and  uses  subroutine 
<RAYINT>  and  subroutine  <APPROX>. 

Output  files; 

RIes  <F_matr>  -  arrays  of  matric  of  corresponding  approximation 
for  each  receiver 

RIe  <Fparam>  -  parameters  of  matrices  for  each  receiver 
RIes  <FJnt>  -  results  of  multiplications  of  calculated  matrix 
and  model  structure 

RIe  <F_st>  -  array  (LST)  of  number  of  all  rays  which  cross  the 
corresponding  discret  of  model 
Compilation 

mf77  mat_nl.for  -ol  -486 
RUN 

ndprun  mat_nl.ltl 


2.3  The  solution  of  the  inverse  probiem  for  satellite  nonlinear 

Radiotomography. 


In  this  section  we  will  cite  the  results  of  solution  of  inverse  nonlinear  RT 
problem.  As  we  noted  earlier,  equation  systems  (9-10)  are  nonlinear  equatons 
relatively  to  electron  density  distribution,  because  electron  density  distribution 
determines  also  the  structure  of  operators  A  and  L  themselves  in  (9-10).  The 
successive  approximation  method  is  the  most  suitable  method  to  solve  such  linear 
equation  systems.  This  approach  is  correct  because  the  ray  distortion  is  not  strong  in 
the  majority  of  practically  interesting  cases  and  operators  of  linear  and  nonlinear  RT 
are  close.  On  the  first  stage  of  successive  approximation  method  matrixes  L  and  A 
are  calculated  using  the  initial  approximation  N.  Then  one  of  the  linear  equation 
systems  is  solving  (9  or  10  correspondly)  and  the  first  approximation  of  electron 
density  distribution  N  is  finding.  On  the  second  stage  the  foliowing  approximations  for 
matrixes  A  or  L  are  calculating  using  the  first  approximation  N.  Then  again  the 
corresponding  linear  equation  system  is  solving  and  the  second  approximation  of 
electron  density  distribution  N  is  finding.  Further  this  process  can  be  repeated.  Thus, 
here,  on  every  stage  ,  the  linear  equation  system  is  solving,  using  the  operator  matrix, 
which  was  constructed  according  to  previous  approximation  N. 

There  ixist  a  lot  of  variants  of  solution  of  linear  RT  problem,  which  give  different 
results.  To  illustrate  the  solution  of  nonlinear  RT  problem  it  is  correct  to  use  the  same 
method  in  a  case  of  linear  RT,  to  make  it  possible  to  distinguish  the  difference 
between  linear  and  nonlinear  RT  for  our  analysis.  In  the  following  examples,  the 
method  C  of  building  the  operators  of  direct  problem  was  used.  (C  is  built  with  the 
linear  product  approximation).  The  use  of  other  methods  of  building  the  operators  of 
direct  problem,  the  same  for  linear  and  nonlinear  case,  gives  the  same  results.  Here, 
to  solve  the  RT  problems,  the  phase  9  and  phase-difference  or  Doppler  (10) 
approach  will  be  used.  Though  the  use  of  phase  approach  over  TEC  Is  appeared  to 
be  unperspective  (  we  repeatedly  criticized  this  method),  however,  here,  it  is  usefull  to 
cite  the  results  of  the  influence  of  refractable  corrections  on  both  approaches. 

To  estimate  the  power  of  the  influence  of  refractable  corrections  on  the 
solution  of  RT  problems  we  will  use  several  simple  ionospheric  models.  It  is  obvious, 
that  the  qualitative  sight  of  the  influence  of  refractable  corrections  on  the  solution  of 
RT  problems  slightly  depends  on  the  model,  and  it  (sight)  is  determined  mainly  by 
the  parameter  aAN  /  (9,10)  which  characterizes  the  variations  of  electron  density. 

As  a  first  model,  the  model  of  quasilayered  ionosphere  with  a  trough  and  two 
irregularities  on  the  border  of  the  trough  (fig.25),  with  maximum  of  electron  density 

equal  to  (f=150  Mhz)  was  used.  Phase  curvatures  of  this  model  were 

already  cited  in  the  previous  section  (fig.  19-34).  Fig. 26-28  show  the  radiotomographic 
reconstruction  results  of  given  model;  fig. 26  -  reconstruction  by  linear  RT  method, 


fig. 27  -  reconstruction  by  nonlinear  RT  method  after  the  first  stage,  fig. 28  - 
reconstruction  by  nonlinear  RT  method  after  the  second  stage.  It  is  notable  that  the 
reconstruction  quality  became  better.  This  is  also  demonstrated  by  reconstruction 
mistakes,  for  example,  for  linear  RT  the  mistakes  {S2  =  12.4%  and 

^00  =  37.9%)^  and  for  nonlinear  RT  after  the  first  stage  {S2  =  8.4%  and 
=  20.8%).  Also  the  congruence  of  the  magnitude  of  electron  density  maximum  for 

nonlinear  case  5.08  (in  units  against  5.59  in  linear  case  is  essentially  better. 

This  essential  improvement  of  reconstructions  is  not  amasing,  because  the  operator 
approximation  of  direct  problem  essentially  improves  at  high  electron  density 
A./V  «  5*10^^  with  the  use  of  nonlinear  RT  approach  (fig. 21 -24).  In  a  case  of  lower 
electron  density  A#  »  0.5*10^^  the  operators  approximation  of  direct  problem  are 
closed  for  linear  and  nonlinear  RT  (fig.  19-20)  and  the  solution  results  of  inverse 
problem  practically  coincide. 

As  a  second  ionospheric  model  the  model  of  quasiwave  structure  on  the 
background  of  the  layered  ionosphere  (fig. 29)  with  electron  density  maximum 
was  used.  On  such  maximum  level  the  results  of  RT  reconstructions  by 
linear  and  nonlinear  RT  methods  are  closed,  however  the  mistakes  in  nonlinear  case 
(<52  =  9.5%  and  <5„  =  24.4%)  are  slightly  less  that  the  linear  case  mistakes 
(62  =  10%  and  =  27.2%).  The  closeness  of  RT  reconstruction  results  is  explained 
by  quite  good  linear  operator  approximation,  what  is  illustrated  in  fig. 32,  where  the 
real  phase  is  illustrated  by  Intire  line,  calculated  with  the  respect  to  retractable 
corrections  and  the  result  of  product  of  operator  (9)  linear  approximation  and  the 
model  is  illustrated  by  dashed  line.  The  situations  in  reconstructions  notably  changes 
with  the  increasing  of  electron  density  maximum  of  quasiwave  structure  to  5*10*^ 

(fig. 33).  The  difference  between  linear  and  nonlinear  RT  now  is  becoming  more 
essential.  The  reconstruction  quality  received  by  means  of  linear  RT  method  is  not 
satisfying:  in  fig.34  it  is  difficult  to  guess  where  the  quasiwave  structures  are 
(<^2  =  13.4%  and  =  32.5%).  The  result  of  the  first  stage  of  nonlinear  RT 
reconstruction  is  cited  in  fig.35,  the  second  stage  in  fig.36.  Seen  is,  with  the  use  of 
nonlinear  RT  approach  the  situation  became  notably  better,  what  is  demonstrated  by 
the  reconstruction  mistakes  {S2  =  8.19%  and  6^  =  24.4%)  after  the  first  stage, 
(<^2  =8%  and  =24%)  after  the  second  stage.  In  the  last  case  the  difference 
between  linear  (  dashed  line  in  fig.37-38)  and  nonlinear  (entire  line,  which  practically 
coincides  with  real  phase)  approximation  is  more  essential  compared  with  the  less 


maximum,  what  is  seen  well  in  fig. 37-38,  where  this  difference  is  cited  for  first  two 
receivers. 

The  results,  cited  earlier,  are  related  to  reconstructions  over  phase  and  TEC. 
As  it  would  be  seen  lower,  the  same  results  in  a  meaning  of  influence  of  refractable 
corrections,  are  obtaining  also  under  the  reconstruction  over  phase  difference  (TEC 
or  Doppler  difference).  In  fig. 39  shown  is  the  follovnng  ionospheric  model  with  a 
trough,  possitive  irregularity  from  the  left  border  of  trough  and  negative  irregularity 
from  the  right .  The  lower  reconstruction  quality  in  linear  approach  according  to  phase 
measurings  (S2  =  12.7%  and  =  38.4%)  is  illustrated  by  fig.40,  where  strong 
artefacts  and  distortions  do  not  allow  to  select  Initial  irregularities.  The  use  of 
sequential  stages  of  nonlinear  RT  method  over  phase  differences  gives  better  results: 

I  stage  (fig.4l,  ^2  =  9.3%  and  =  46.5%),  II  stage  (fig.42, 
S2  -  9.18%  and  =  48.8%),  III  stage  (fig.43,  S2  =  9.1%  and  =  46.7%),  IV  stage 
(fig.44,  S2  =  9.12%  and  =  46.8%),.  The  use  of  three  stages  is  enough  here 

as  the  quality  improvement  doesn’t  take  place  with  the  following  stages.  The  use  of 
Doppler  method  essentially  improves  the  reconstruction  results  both  in  linear  and 
nonlinear  cases,  what  is  explained  by  higher  sensitivity  [25]  of  difference  Doppler 
method.  The  Doppler  reconstruction  results  are  given  in  fig.45  (linear  case, 
82  =  7%  and  5^  =  20%)  and  fig.46  (nonlinear  case,  82  =  5.6%  and  8^  =  19%)  , 
where  high  reconstruction  quality  is  seen  well. 

The  last  fourth  model  which  illustrates  the  nonlinear  RT  methods  is  shown  in 
fig.47.  It  is  a  quite  complete  structure  with  a  trough  and  four  irregularities  (  two 
"positive"  r  « 1300^  and  two  "negative"  r  «-50/:OT,r  » lOOOA/w).  Such 

combinations  of  "positive"  and  "negative"  irregularities,  overlaping  each  other  for 
different  receivers,  usually  are  the  most  difficult  for  RT  reconstructions.  The  electron 
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density  level  is  supposed  to  be  high  enough  AN  5*10  m  ,  to  make 
refractable  effects  exhibit  notably.  The  result  of  RT  reconstruction  by  the  linear 
Doppler  RT  method  is  shown  in  fig.48.  (S2  =  7.8%  and  8^0  —  18.7%).  Here,  the 
possitive  irregularities  are  seen  well,  but  the  negative  practically  are  not  exhibited.  The 
results  of  RT  reconstruction  made  by  nonlinear  Doppler  RT  method  are  much  better, 
the  are  shown  in  fig.40  (  after  the  first  stage  82  =  6.5%  and  8^  =  16.5%)  and  in 
fig.41  (  after  the  second  stage  82  =  5.5%  and  8^  =  13.8%).  The  attention  should 
be  payed  to  quite  high  reconstruction  quality  in  fig.41,  where  the  main  qualitative 
pecularities  of  reconstructed  structure  are  exhibited  well.  Let  us  note  that  here  and 
earlier  during  the  step  by  step  RT  reconstruction  only  small  quantity  of  iterations  2-4 
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by  ART  and  MART  methods  were  used.  The  mistake  ,  compared  with  the  mistake  of 
operator  approximation  souid  be  ieft  in  the  first  part  of  equations  (9-10).  Than,  if  we 
increase  the  iteration  quantity  and  dicrease  the  mistake  in  the  right  part  ( the  mistake 
according  to  data  I  or  D  )  the  solution  will  go  away  from  the  real  and  they  (solutions) 
will  diverge.  Than  the  matrix  operator  calculation  for  the  next  stage  will  have  the  large 
mistakes  and  the  sequential  stage  will  not  lead  to  the  improvement  of  results. 
However  such  a  small  iteration  quantity  gives  usually  good  coincidence  of  the  right 
part  of  equations  (9-10).  As  an  example,  let  us  cite  the  result  of  comparison  of 
Doppler  frequencies  for  the  second  and  the  third  receivers  after  only  the  two  iterations 
of  the  second  stage  (  relative  mistakes  of  Doppler  frequencies 
S2  =  11%  and  S^o  =  24%).  Fig. 51  and  fig.52  show  the  real  Doppler  frequencies 
for  by  entire  line,  and  Doppler  frequencies  after  RT  reconstructions  from  fig. 50  -  by 
dashed  line. 

Generally  it  should  be  stated  that  the  retractable  effects  becomes  essential 
when  the  deflection  of  retractable  rays  from  direct  becomes  comparable  with  the  size 
of  element  digitization  in  RT  problem.  The  deflections  of  retractable  rays  are 
determined  mainly  by  electron  density  difference  AN.  On  the  probing  frequencies 

equal  to  hundreds  of  MHz,  usually,  the  refraction  is  ignored  for  AN  <  1*10*^ and 
for  Size  of  element  digitization  more  than  20  km.  The  refraction  should  be  respected 

for  large  AN  >(2-4)*10’^/w“^,  for  example,  in  the  years  of  solar  activity,  and  also  for 
probing  on  frequencies  equal  to  50-100  MHz.  More  precise  estimations  ,  for 
intermediate  cases,  can  be  obtained  from  6-8,  and  also  by  computer  modelling  with 
the  help  of  programs  presented  here. 

Description  of  the  program  for  solution  of  inverse  problem  Nonlinear  RT 

Program  <slv_nl.for> 

This  program  solves  inverse  problem  for  TEC,  Phase  or  Phase-difference  RT 
by  means  of  different  algorithms:  ART,  MART,  SIRT.  It  calculates  errors  of 
reconstruction  in  metric  C  and  L2  in  dependence  upon  data  errors  and  noises. 

Input  parameters  and  files. 

When  staring  program,  you  may  choose  the  type  of  RT  measurements 
Answer  =  1  -  Phase  RT:  Tec_ph  =  1  — >  TEC  measurements 
Tec_ph  =  2  — >  Phase  measurements 
Answer  =  2  :  Phase-difference  RT 
Answer,  Tec_ph  are  input  from  the  screen 


NF,  NR  -  number  of  discrets  on  the  horizontal  and  vertical  grid 
NREC  -  number  of  receivers 

T_MOD'  -  input  file  of  model  structure  (program  <dir_rth1.for>) 

'F_XO’  -  input  file  of  initial  guess 

(program  <dir_rth1.for>  or  function  <guess>) 

'Fparam'  -  input  file  with  parameters  of  matrix 
XIST  -  model  structure  from  input  file  'F_MOD' 

File  'name_F.slv'  containes  Nrec  lines  with  names  of  input  and  output  files. 

The  example  of  this  file  is  enclosed. 

line  1  -  name  of  file  with  parameters  of  matrix  (Fparam) 

line  2:  name!  -  name  of  input  file  with  matrix 

name2  -  name  of  input  file  with  either  TEC  or  TEC-difference 
('Finput') 

names  -  name  of  output  file  with  results  of  multiplications 
of  input  matric  and  reconstructed  structure  ('Fout') 
line  Nrec+2:  name  of  file  with  number  of  all  rays  which  cross 
the  corresponding  discret  of  model 
line  Nrec+3:  name  file  with  model,  name  file  with  guess, 
name  file  with  reconstruction 
line  Nrec+4:  name  file  with  errors 
Npi  -  array  of  contants  (2pn)  for  each  receiver 
er_n  -  level  of  noise  in  % 

Npi,  er_n  are  input  from  the  screen 

Nsolve  -  the  method  of  solution,  it  is  input  from  the  screen 

REL  -  array  of  relaxation  parameter 

LMAX  -  max  number  of  nonzero  matric  elements  for  each  receiver 
AZ  -  array  of  nonzero  matric  elements  for  all  receivers 
NST  -  array  of  corresponding  column  nonzero  matric  elements  for 
all  receivers 

LST  -  array  of  number  of  all  rays  which  cross  the 

corresponding  discret  of  model  (SIRT  algorithm)  from 
input  file  <F_ST>,  it  is  similar  to  LST  in  program  <mat_nl.for> 

Niter  -  number  of  iterations,  it  is  input  from  screen 
Subroutine  <DEFSYS>  -  solution  of  linear  system  equations  for 
one  ray. 

Subroutine  <ercl2>  calculates  errors  in  metric  C  and  L2 
Subroutine  <VMINM>  calculates  min  and  max  of  array 


Function  <RAN>  calculates  random  values  in  [0,1] 

Function  <GUESS>  -  for  initial  guess 

RMAX,  RM,  Zmax,  ZSM,  B1,  B2  -  parameters  for  function  <GUESS> 
Output  files: 

File  <F_REC>  -  file  with  the  reconstruction 
File  <er_solv>  -  errors  of  right  items  and  reconstruction  in 
metric  C  and  L2 

Files  <Fout>  -  results  of  multiplications  of  the  metric 
and  reconstructed  structure 

Execution 

mf77  slv_nl.for  -ol  -486 
RUN 

ndprun  slv_nl.ltl 
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